sandbox/symposium2013/optimize.HFindexes.R
In braverock/PortfolioAnalytics: Portfolio Analysis, Including Numerical Methods for Optimization of Portfolios
### For Presentation at FactSet's 2013 US Investment Process Symposium
# November 10 - 12 , 2013
# Peter Carl
### Load the necessary packages
# Include optimizer packages
require(PortfolioAnalytics)
require(quantmod)
library(DEoptim)
library(ROI)
require(ROI.plugin.quadprog)
require(ROI.plugin.glpk)
# ... and multi-core packages
require(foreach)
require(doMC)
registerDoMC(6)
# Available on r-forge
# require(FactorAnalytics) # development version > build
### Set script constants
runname='historical.moments'
rebalance_period = 'years' #'quarters' # uses endpoints identifiers from xts; how to do semi-annual?
clean = "boudt" # "none"
permutations = 2000
p=1-1/12 # set confidence for VaR/mETL for monthly data
### Description
# Seven assets, in this case hedge fund indexes representing different styles and two 'types':
# Convertible Arbitrage | Non-directional
# Equity Market Neutral | Non-directional
# Fixed Income Arbitrage | Non-directional
# Event Driven | Non-directional
# CTA Global | Directional
# Global Macro | Directional
# Long/Short Equity | Directional
# see analyze.HFindexes.R for more detail
# Set up seven objectives as buoy portfolios:
# - Equal contribution to...
# 1 Weight
# 2 Variance
# 3 Risk (mETL)
# - Reward to Risk ratio of...
# 4 Mean-Variance
# 5 Mean-mETL
# - Minimum...
# 6 Variance
# 7 Risk (mETL)
# Add constraints
# - Box constraints - 5% to 30%
# - Group constraints - Non-directional constrained to 20-70%; Directional between 10-50%
# - Rebalancing period - quarterly
# - Turnover constraints - TBD
#------------------------------------------------------------------------
# Set up an initial portfolio object with constraints and objectives using
# v2 specification
# Create initial portfolio object used to initialize ALL the bouy portfolios
init.portf <- portfolio.spec(assets=colnames(R),
weight_seq=generatesequence(by=0.005)
)
# Add leverage constraint
init.portf <- add.constraint(portfolio=init.portf,
type="leverage"
)
# Add box constraint
init.portf <- add.constraint(portfolio=init.portf,
type="box",
min=0.05,
max=0.3
)
# Add group constraint
init.portf <- add.constraint(portfolio=init.portf, type="group",
groups=list(c(1:4),
c(5:7)),
group_min=c(0.25,.05),
group_max=c(0.85,0.75)
)
# print(init.portf)
# summary(init.portf)
#------------------------------------------------------------------------
### Construct BUOY 1: Constrained Mean-StdDev Portfolio - using ROI
# Add the return and sd objectives to the constraints created above
MeanSD.portf <- add.objective(portfolio=init.portf,
type="return", # the kind of objective this is
name="mean" # name of the function
)
MeanSD.portf <- add.objective(portfolio=MeanSD.portf,
type="risk", # the kind of objective this is
name="var", # name of the function
arguments=list(clean=clean)
)
### Construct BUOY 2: Constrained Mean-mETL Portfolio - using ROI
#@ Cannot maximize mean return per unit ETL with ROI, consider using
#@ random portfolios or DEoptim. - RB
# Add the return and mETL objectives
MeanmETL.portf <- add.objective(portfolio=init.portf,
type="return", # the kind of objective this is
name="mean" # name of the function
, multiplier=-12
)
MeanmETL.portf <- add.objective(portfolio=MeanmETL.portf,
type="risk", # the kind of objective this is
name="ES", # the function to minimize
arguments=list(p=p, clean=clean)
)
### Construct BUOY 3: Constrained Minimum Variance Portfolio - using ROI
# Add the variance objective
MinSD.portf <- add.objective(portfolio=init.portf,
type="risk", # the kind of objective this is
name="var", # name of the function
arguments=list(p=p, clean=clean)
)
### Construct BUOY 4: Constrained Minimum mETL Portfolio - using ROI
# Add the mETL objective
MinmETL.portf <- add.objective(portfolio=init.portf,
type="risk", # the kind of objective this is
name="ES", # the function to minimize
arguments=list(p=p, clean=clean)
)
### Construct BUOY 5: Constrained Equal Variance Contribution Portfolio - using RP
#@ - Add the sub-objectives first. Adding these 3 objectives means that we are
#@ maximizing mean per unit StdDev with equal volatility contribution portfolios. - RB
# Without a sub-objective, we get a somewhat undefined result, since there are (potentially) many Equal SD contribution portfolios.
# MRCSD.portf <- add.objective(portfolio=init.portf,
# type="risk",
# name="StdDev"
# ) # OR
# MRCSD.portf <- add.objective(portfolio=init.portf,
# type="return",
# name="mean"
# )
MRCSD.portf <- add.objective(portfolio=init.portf,
type="risk_budget",
name="StdDev",
min_concentration=TRUE,
arguments = list(clean=clean)
)
# MRCSD.portf$constraints[[1]]$min_sum = 0.99 # set to speed up RP
# MRCSD.portf$constraints[[1]]$max_sum = 1.01
### Construct BUOY 6: Constrained Equal mETL Contribution Portfolio - using RP
#@ Add the sub-objectives first. These should be added to the MRCmETL portfolio.
#@ All objectives below mean that we are maximizing mean return per unit ES with
#@ equal ES contribution. - RB
# Without a sub-objective, we get a somewhat undefined result, since there are (potentially) many Equal SD contribution portfolios.
# MRCmETL.portf <- add.objective(portfolio=init.portf,
# type="risk",
# name="ES"
# ) # OR
MRCmETL.portf <- add.objective(portfolio=init.portf,
type="return",
name="mean"
)
MRCmETL.portf <- add.objective(MRCmETL.portf,
type="risk_budget",
name="ES",
min_concentration=TRUE,
arguments = list(p=p, clean=clean)
)
# Calculate portfolio variance, but don't use it in the objective; used only for plots
MRCmETL.portf <- add.objective(portfolio=MRCmETL.portf,
type="risk", # the kind of objective this is
name="StdDev", # the function to minimize
enabled=TRUE, # enable or disable the objective
multiplier=0, # calculate it but don't use it in the objective
arguments=list(clean=clean)
)
# MRCmETL.portf$constraints[[1]]$min_sum = 0.99 # set to speed up RP
# MRCmETL.portf$constraints[[1]]$max_sum = 1.01
### Construct BUOY 7: Equal Weight Portfolio
# There's only one, so create a portfolio object with all the objectives we want calculated.
EqWgt.portf <- portfolio.spec(assets=colnames(R))
EqWgt.portf <- add.constraint(portfolio=EqWgt.portf, type="leverage", min_sum=1, max_sum=1)
EqWgt.portf <- add.objective(portfolio=EqWgt.portf, type="return", name="mean")
EqWgt.portf <- add.objective(portfolio=EqWgt.portf, type="risk_budget", name="ES", arguments=list(p=p, clean=clean))
EqWgt.portf <- add.objective(portfolio=EqWgt.portf, type="risk_budget", name="StdDev", arguments=list(clean=clean))
### Construct BUOY 8: Inverse Volatility Portfolio
# There's only one, so create a portfolio object with all the objectives we want calculated.
VolWgt.portf <- portfolio.spec(assets=colnames(R))
VolWgt.portf <- add.constraint(portfolio=VolWgt.portf, type="leverage", min_sum=0.99, max_sum=1.01)
VolWgt.portf <- add.objective(portfolio=VolWgt.portf, type="return", name="mean")
VolWgt.portf <- add.objective(portfolio=VolWgt.portf, type="risk_budget", name="ES", arguments=list(p=p, clean=clean))
VolWgt.portf <- add.objective(portfolio=VolWgt.portf, type="risk_budget", name="StdDev", arguments=list(clean=clean))
# REMOVED - to much to show already
# ### Construct RISK BUDGET Portfolio
# ConstrConcmETL.portf <- portfolio.spec(assets=colnames(R),
# weight_seq=generatesequence(by=0.005)
# )
# # Add leverage constraint
# ConstrConcmETL.portf <- add.constraint(portfolio=RiskBudget.portf,
# type="leverage",
# min_sum=0.99, # set to speed up RP, DE
# max_sum=1.01
# )
# # Establish position bounds
# ConstrConcmETL.portf <- add.constraint(portfolio=ConstrConcmETL.portf,
# type="box",
# min=0.01, # leave relatively unconstrained
# max=1.0
# )
# # Maximize mean return
# ConstrConcmETL.portf <- add.objective(portfolio=ConstrConcmETL.portf,
# type="return", # maximize return
# name="mean",
# multiplier=12
# )
# # Add a risk measure
# # Use ES to be consistent with risk measures in other BUOY portfolios
# ConstrConcmETL.portf <- add.objective(portfolio=ConstrConcmETL.portf,
# type="risk",
# name="ETL", # using a different name to avoid clobbering slot below, workaround for bug
# multiplier=1,
# arguments = list(p=p, clean=clean)
# )
#
# # Set contribution limits
# ConstrConcmETL.portf <- add.objective(portfolio=ConstrConcmETL.portf,
# type="risk_budget",
# name="ES",
# max_prisk=0.3, # Sets the maximum percentage contribution to risk
# arguments = list(p=p, clean=clean)
# )
# Calculate portfolio variance, but don't use it in the objective; used only for plots
# ConstrConcmETL.portf <- add.objective(portfolio=ConstrConcmETL.portf,
# type="risk", # the kind of objective this is
# name="StdDev", # the function to minimize
# enabled=TRUE, # enable or disable the objective
# multiplier=0, # calculate it but don't use it in the objective
# arguments=list(clean=clean)
# )
#------------------------------------------------------------------------
### Evaluate portfolio objective objects
#------------------------------------------------------------------------
# Generate a single set of random portfolios to evaluate against all RP constraint sets
print(paste('constructing random portfolios at',Sys.time()))
# Modify the init.portf specification to get RP running
rp.portf <- init.portf
# rp.portf$constraints[[1]]$min_sum = 0.99 # set to speed up RP
# rp.portf$constraints[[1]]$max_sum = 1.01
rp.portf$constraints[[1]]$min_sum = 1.00 # for more accuracy
rp.portf$constraints[[1]]$max_sum = 1.00
# rp = random_portfolios(portfolio=rp.portf, permutations=30000, max_permutations=400) # will get fewer with less accuracy
load(file=paste(resultsdir,'random-portfolios-2013-10-05.historical.moments.rda'))
rp.mean = apply(rp, 1, function(w) mean(R %*% w))
rp.sd = apply(rp, 1, function(x) StdDev(R=R, weights=x, p=p, clean=clean))
plot(rp.sd, rp.mean, col="darkgray", cex=0.5)
# REMOVED: This was fruitless
# rp1 = random_portfolios(portfolio=rp.portf, permutations=10000, max_permutations=400, rp_method="sample")
# rp1.mean = apply(rp1, 1, function(w) mean(R %*% w))
# rp1.sd = apply(rp1, 1, function(x) StdDev(R=R, weights=x, p=p))
# rp1.etl=NULL; for(i in 1:NROW(rp1)) {rp1.etl[i]=ETL(R=R, weights=as.vector(rp1[i,]), p=p, portfolio_method="component")[[1]]}
# plot(rp1.sd, rp1.mean, col="gray", cex=0.5)
#
# rp2 = random_portfolios(portfolio=rp.portf, permutations=10000, max_permutations=400, rp_method="simplex", fev=2)
# rp2.mean = apply(rp2, 1, function(w) mean(R %*% w))
# rp2.sd = apply(rp2, 1, function(x) StdDev(R=R, weights=x, p=p))
# points(rp2.sd,rp2.mean, col="blue", cex=0.5)
#
# rp3 = random_portfolios(portfolio=rp.portf, permutations=10000, max_permutations=400, rp_method="grid")
# rp3.mean = apply(rp3, 1, function(w) mean(R %*% w))
# rp3.sd = apply(rp3, 1, function(x) StdDev(R=R, weights=x, p=p))
# points(rp3.sd,rp3.mean, col="green", cex=0.5)
# print(paste('done constructing random portfolios at',Sys.time()))
# save(rp,file=paste(resultsdir, 'random-portfolios-', Sys.Date(), '-', runname, '.rda',sep=''))
start_time<-Sys.time()
print(paste('Starting optimization at',Sys.time()))
### Evaluate BUOY 1: Constrained Mean-StdDev Portfolio - with ROI
MeanSD.ROI<-optimize.portfolio(R=R,
portfolio=MeanSD.portf,
optimize_method='ROI',
trace=TRUE
)
plot(MeanSD.ROI, risk.col="StdDev", return.col="mean", rp=permutations, chart.assets=TRUE, main="Mean-Volatility Portfolio")
save(MeanSD.ROI,file=paste(resultsdir, 'MeanSD-', Sys.Date(), '-', runname, '.rda',sep='')) # Save the results
print(paste('Completed meanSD optimization at',Sys.time(),'moving on to meanmETL'))
### Evaluate BUOY 2: Constrained Mean-mETL Portfolio - with DE or RND
# Not possible with ROI - RB
# MeanmETL.ROI<-optimize.portfolio(R=R,
# portfolio=MeanmETL.portf,
# optimize_method='ROI',
# trace=TRUE, verbose=TRUE
# )
#
# So use random portfolios instead
MeanmETL.RND<-optimize.portfolio(R=R,
portfolio=MeanmETL.portf,
optimize_method='random',
rp=rp,
trace=TRUE
)
plot(MeanmETL.RND, risk.col="StdDev", return.col="mean", rp=permutations, chart.assets=TRUE, main="Mean-mETL Portfolio")
plot(MeanmETL.RND, risk.col="ES", return.col="mean", rp=permutations, chart.assets=TRUE, main="Mean-mETL Portfolio")
save(MeanmETL.RND,file=paste(resultsdir, 'MeanETL-RP-', Sys.Date(), '-', runname, '.rda',sep=''))
chart.EfficientFrontier(MeanmETL.RND)
# OR with DE optim
MeanmETL.DE<-optimize.portfolio(R=R,
portfolio=MeanmETL.portf,
optimize_method='DEoptim',
search_size=20000,
initialpop=rp[1:50,] # seed with a starting population that we know fits the constraint space
)
plot(MeanmETL.DE, risk.col="StdDev", return.col="mean", chart.assets=TRUE, main="Mean-mETL Portfolio")
plot(MeanmETL.DE, risk.col="ES", return.col="mean", chart.assets=TRUE, main="Mean-mETL Portfolio")
save(MeanmETL.DE,file=paste(resultsdir, 'MeanETL-DE-', Sys.Date(), '-', runname, '.rda',sep=''))
print(paste('Completed meanmETL optimization at',Sys.time(),'moving on to MinSD'))
### Evaluate BUOY 3: Constrained Minimum Variance Portfolio - with ROI
MinSD.ROI<-optimize.portfolio(R=R,
portfolio=MinSD.portf,
optimize_method='ROI',
trace=TRUE
) #
plot(MinSD.ROI, risk.col="StdDev", return.col="mean", rp=permutations, chart.assets=TRUE, main="Minimum Volatility Portfolio with ROI")
save(MinSD.ROI,file=paste(resultsdir, 'MinSD-', Sys.Date(), '-', runname, '.rda',sep=''))
print(paste('Completed MinSD optimization at',Sys.time(),'moving on to MinmETL'))
# OR with random portfolios
# MinSD.RND<-optimize.portfolio(R=R,
# portfolio=MinSD.portf,
# optimize_method='random',
# rp=rp,
# trace=TRUE, verbose=TRUE
# ) #
# plot(MinSD.RND, risk.col="StdDev", return.col="mean", chart.assets=TRUE, main="Minimum Volatility Portfolio with RP")
### Evaluate BUOY 4: Constrained Minimum mETL Portfolio - with ROI
MinmETL.ROI<-optimize.portfolio(R=R,
portfolio=MinmETL.portf,
optimize_method='ROI',
trace=TRUE, verbose=TRUE
)
plot(MinmETL.ROI, risk.col="StdDev", return.col="mean", rp=permutations, chart.assets=TRUE, main="Minimum mETL Portfolio")
plot(MinmETL.ROI, risk.col="ES", return.col="mean", rp=permutations, chart.assets=TRUE, main="Minimum mETL Portfolio")
save(MinmETL.ROI,file=paste(resultsdir, 'MinmETL-', Sys.Date(), '-', runname, '.rda',sep=''))
print(paste('Completed MinmETL optimization at',Sys.time(),'moving on to MRCSD'))
### Evaluate BUOY 5: Constrained Equal Variance Contribution Portfolio - with RP
# MRCSD.RND<-optimize.portfolio(R=R,
# portfolio=MRCSD.portf,
# optimize_method='random',
# rp=rp,
# trace=TRUE
# )
# plot(MRCSD.RND, risk.col="StdDev", return.col="mean", chart.assets=TRUE, main="Equal Volatility Contribution Portfolio")
# chart.RiskBudget(MRCSD.RND, risk.type="percentage", neighbors=25)
# save(MRCSD.RND,file=paste(resultsdir, 'MRCSD.RND-', Sys.Date(), '-', runname, '.rda',sep=''))
# ... not a very satisfying solution
# OR DE optim - this gets very close (a nice, straight line), so use it
MRCSD.DE<-optimize.portfolio(R=R,
portfolio=MRCSD.portf,
optimize_method='DEoptim',
search_size=20000,
itermax=400,
initialpop=rp[1:50,], # seed with a starting population that we know fits the constraint space
trace=FALSE
)
plot(MRCSD.DE, risk.col="StdDev", return.col="mean", chart.assets=TRUE, main="Equal Volatility Contribution Portfolio")
chart.RiskBudget(MRCSD.DE, risk.type="percentage", neighbors=25)
save(MRCSD.DE,file=paste(resultsdir, 'MRCSD.DE-', Sys.Date(), '-', runname, '.rda',sep=''))
print(paste('Completed MRCSD optimization at',Sys.time(),'moving on to MRCmETL'))
### Evaluate BUOY 6: Constrained Equal mETL Contribution Portfolio - with RP
MRCmETL.RND<-optimize.portfolio(R=R,
portfolio=MRCmETL.portf,
optimize_method='random',
rp=rp,
trace=TRUE
) #
plot(MRCmETL.RND, risk.col="StdDev", return.col="mean", chart.assets=TRUE, main="Equal mETL Contribution Portfolio")
plot(MRCmETL.RND, risk.col="ES", return.col="mean", chart.assets=TRUE, main="Equal mETL Contribution Portfolio")
chart.RiskBudget(MRCmETL.RND, neighbors=25)
save(MRCmETL.RND,file=paste(resultsdir, 'MRCmETL-', Sys.Date(), '-', runname, '.rda',sep=''))
# OR DE optim -
MRCmETL.DE<-optimize.portfolio(R=R,
portfolio=MRCmETL.portf,
optimize_method='DEoptim',
search_size=20000,
NP=200,
initialpop=rp[1:50,], # seed with a starting population that we know fits the constraint space
trace=FALSE
)
plot(MRCmETL.DE, risk.col="StdDev", return.col="mean", chart.assets=TRUE, main="Equal Volatility Contribution Portfolio")
chart.RiskBudget(MRCmETL.DE, risk.type="percentage", neighbors=25)
save(MRCmETL.DE,file=paste(resultsdir, 'MRCmETL.DE-', Sys.Date(), '-', runname, '.rda',sep=''))
# # test it unconstrained:
# unconstr.portf <- portfolio.spec(assets=colnames(R),
# weight_seq=generatesequence(by=0.005)
# )
# unconstr.portf <- add.constraint(portfolio=unconstr.portf,
# type="leverage",
# min_sum=0.99, # set to speed up RP
# max_sum=1.01
# )
# # Establish position bounds
# unconstr.portf <- add.constraint(portfolio=unconstr.portf,
# type="box",
# min=0.01,
# max=1.0
# )
# MRCmETLun.portf <- add.objective(portfolio=unconstr.portf,
# type="return",
# name="mean"
# )
# MRCmETLun.portf <- add.objective(MRCmETL.portf,
# type="risk_budget",
# name="ES",
# min_concentration=TRUE,
# arguments = list(p=p, clean=clean)
# )
#
# # ...in DE optim -
# MRCmETLun.DE<-optimize.portfolio(R=R,
# portfolio=MRCmETLun.portf,
# optimize_method='DEoptim',
# search_size=20000,
# NP=200,
# initialpop=rp[1:50,], # seed with a starting population that we know fits the constraint space
# trace=FALSE
# )
print(paste('Completed MRCmETL optimization at',Sys.time(),'moving on to RiskBudget'))
### Evaluate BUOY 7: Equal Weight Portfolio
# Calculate the objective measures for the equal weight portfolio
EqWgt.opt <- equal.weight(R=R, portfolio=EqWgt.portf)
### Evaluate BUOY 8: Inverse Volatility Portfolio
volatility.weight <- function (R, portfolio, ...)
{
if (!is.portfolio(portfolio))
stop("portfolio object passed in must be of class 'portfolio'")
assets <- portfolio$assets
nassets <- length(assets)
weights <- as.vector((1/StdDev(R))/sum(1/StdDev(R)))
names(weights) <- names(assets)
if (ncol(R) != nassets) {
if (ncol(R) > nassets) {
R <- R[, 1:nassets]
warning("number of assets is less than number of columns in returns object, subsetting returns object.")
}
else {
stop("number of assets is greater than number of columns in returns object")
}
}
out <- constrained_objective(w = weights, R = R, portfolio = portfolio,
trace = TRUE, ...)$objective_measures
return(structure(list(R = R, weights = weights, objective_measures = out,
call = match.call(), portfolio = portfolio), class = c("optimize.portfolio.invol",
"optimize.portfolio")))
}
# Calculate the objective measures for the vol weight portfolio
VolWgt.opt <- volatility.weight(R=R, portfolio=VolWgt.portf)
# REMOVED
# ### Evaluate Constrained Concentration to mETL Portfolio - with DE
# # registerDoSEQ() # turn off parallelization to keep the trace data
# ConstrConcmETL.DE<-optimize.portfolio(R=R,
# portfolio=ConstrConcmETL.portf,
# optimize_method='DEoptim',
# search_size=40000,
# NP=4000,
# itermax=400,
# trace=FALSE
# )
# # list(c=0.25, # speed of crossover adaption (0,1]
# # CR=0.75) # crossover probability [0,1]
# plot(ConstrConcmETL.DE, risk.col="StdDev", return.col="mean")
# plot(ConstrConcmETL.DE, risk.col="ES", return.col="mean") # several outlier portfolios
# chart.RiskBudget(ConstrConcmETL.DE)
# chart.RiskBudget(ConstrConcmETL.DE, risk.type="percentage")
#
# save(ConstrConcmETL.DE,file=paste(resultsdir, 'ConstrConcmETL-', Sys.Date(), '-', runname, '.rda',sep=''))
# print(ConstrConcmETL.DE$elapsed_time)
print('Done with optimizations.')
#------------------------------------------------------------------------
### Extract data from optimizations for analysis
#------------------------------------------------------------------------
# Combine optimization objects
buoys <- combine.optimizations(list(MeanSD=MeanSD.ROI, MeanmETL=MeanmETL.RND, MinSD=MinSD.ROI, MinmETL=MinmETL.ROI, MRCSD=MRCSD.DE, MRCmETL=MRCmETL.DE, VolWgt=VolWgt.opt, EqWgt=EqWgt.opt))
# chart.Weights(buoys, plot.type="bar", ylim=c(0,1))
#
# #@ Chart the portfolios that have mean and ES as objective measures. - RB
# chart.RiskReward(buoys, risk.col="ES")
# #@ Chart the portfolios that have mean and StdDev as objective measures. - RB
# chart.RiskReward(buoys, risk.col="StdDev")
#
# #@ The MRCmETL and RB optimizations would be good to compare because they are
# #@ similar in that they both include component ES as an objective. - RB
# buoyETL <- combine.optimizations(list(MRCmETL=MRCmETL.RND, RB=RiskBudget.DE, EqWgt=EqWgt.opt))
# chart.RiskBudget(buoyETL, match.col="ES", risk.type="percentage", legend.loc="topright")
#
# #@ Compare the equal weight portfolio and the equal SD contribution portfolio. - RB
# buoyStdDev <- combine.optimizations(list(MRCSD=MRCSD.RND, EqWgt=EqWgt.opt))
# chart.RiskBudget(buoyStdDev, match.col="StdDev", risk.type="absolute", legend.loc="topleft")
Wgts = extractWeights(buoys)
### Extract portfolio measures from each objective
# We can't just extract them, because they aren't all calculated
# so fill them in...
buoys.portfmeas = buoys.contrib.sd = buoys.contrib.es = buoys.perc.sd = buoys.perc.es = NULL
for(i in 1:NROW(Wgts)){
mean = sum(colMeans(R)*Wgts[i,])
sd = StdDev(R, weights=Wgts[i,], portfolio_method="component", clean=clean)
es = ES(R, weights=Wgts[i,], method="modified", portfolio_method="component", p=p, clean=clean)
buoys.portfmeas=rbind(buoys.portfmeas, c(mean, sd[[1]][1], es[[1]][1]))
buoys.contrib.sd= rbind(buoys.contrib.sd,sd[[2]])
buoys.contrib.es= rbind(buoys.contrib.es,es[[2]])
buoys.perc.sd = rbind(buoys.perc.sd,sd[[3]])
buoys.perc.es = rbind(buoys.perc.es,es[[3]])
}
colnames(buoys.portfmeas)=c("Mean", "StdDev", "mETL")
rownames(buoys.portfmeas)=
rownames(buoys.contrib.sd)=
rownames(buoys.contrib.es)=
rownames(buoys.perc.sd) =
rownames(buoys.perc.es) = rownames(Wgts)
colnames(buoys.contrib.sd)=
colnames(buoys.contrib.es)=
colnames(buoys.perc.sd) =
colnames(buoys.perc.es) = colnames(Wgts)
# get the RP portfolios with risk and return pre-calculated
xtract = extractStats(MRCmETL.RND)
save(xtract,file=paste(resultsdir, 'xtract-RPs-', Sys.Date(), '-', runname, '.rda',sep=''))
# columnnames = colnames(xtract)
results.names=rownames(buoys.portfmeas)
# by Asset metrics
assets.portfmeas=as.matrix(scatterFUN(R, FUN="mean"))
assets.portfmeas=cbind(assets.portfmeas, scatterFUN(R, FUN="StdDev", clean=clean))
assets.portfmeas=cbind(assets.portfmeas, scatterFUN(R, FUN="ES", clean=clean))
colnames(assets.portfmeas)=c("Mean", "StdDev", "mETL")
rownames(assets.portfmeas)=colnames(Wgts)
#------------------------------------------------------------------------
# Run select buoy optimizations through time
#------------------------------------------------------------------------
#
# Equal Weight
dates=index(R[endpoints(R, on=rebalance_period)])
EqWgt.w = xts(matrix(rep(1/NCOL(R),length(dates)*NCOL(R)), ncol=NCOL(R)), order.by=dates)
colnames(EqWgt.w)= colnames(R)
EqWgt.R=Return.rebalancing(R, EqWgt.w)
chart.StackedBar(EqWgt.w, colorset=wb13color, gap=0)
VolWgt.w = NULL
for(i in 3:length(dates)){
x = volatility.weight(R=R[paste0("::",dates[i]),], portfolio=VolWgt.portf)
VolWgt.w = rbind(VolWgt.w, x$weights)
}
VolWgt.w = as.xts(VolWgt.w, order.by=dates[-1:-2])
VolWgt.R=Return.rebalancing(R, VolWgt.w)
# Equal SD
MRCSD.DE.t = optimize.portfolio.rebalancing(R=R,
portfolio=MRCSD.portf,
optimize_method='DEoptim',
search_size=2000,
NP=200,
initialpop=rp[1:50,], # seed with a starting population that we know fits the constraint space
trace=FALSE,
rebalance_on=rebalance_period, # uses xts 'endpoints'
trailing_periods=NULL, # calculates from inception
training_period=36) # starts 3 years in to the data history
MRCSD.w = extractWeights(MRCSD.DE.t)
MRCSD.gw = extractGroups(MRCSD.DE.t)
save(MRCSD.DE.t,file=paste(resultsdir, 'MRCSD.DE.t-', Sys.Date(), '-', runname, '.rda',sep=''))
chart.UnStackedBar(MRCSD.w, rotate="horizontal", colorset=wb13color, space=0, las=2)
MRCSD.R=Return.rebalancing(edhec.R, MRCSD.w)
colnames(MRCSD) = "MRCSD"
# Extract perc contrib of mES from results object
x=NULL
for(i in 1:length(names(MRCSD.DE.t))) {
x = rbind(x,MRCSD.DE.t[[i]][["objective_measures"]]$StdDev$pct_contrib_StdDev)
}
MRCSD.DE.pct_contrib_StdDev.t = as.xts(x, order.by=as.POSIXct(names(MRCSD.DE.t)))
chart.UnStackedBar(x.xts, rotate="horizontal", colorset=wb13color, space=0, las=2)
# MRC mETL
MRCmETL.DE.t = optimize.portfolio.rebalancing(R=R,
portfolio=MRCmETL.portf,
optimize_method='DEoptim',
search_size=20000,
NP=200,
initialpop=rp[1:50,], # seed with a starting population that we know fits the constraint space
trace=FALSE,
rebalance_on=rebalance_period, # uses xts 'endpoints'
trailing_periods=NULL, # calculates from inception
training_period=36) # starts 3 years in to the data history
MRCmETL.w = extractWeights(MRCmETL.DE.t)
MRCmETL.gw = extractGroups(MRCmETL.DE.t)
save(MRCmETL.DE.t,file=paste(resultsdir, 'MRCmETL.DE.t-', Sys.Date(), '-', runname, '.rda',sep=''))
chart.UnStackedBar(MRCmETL.w, rotate="horizontal", colorset=wb13color, space=0, las=2)
MRCmETL=Return.rebalancing(edhec.R, MRCmETL.w)
colnames(MRCmETL) = "MRCmETL"
MRCmETL.RND.t = optimize.portfolio.rebalancing(R=R,
portfolio=MRCmETL.portf,
optimize_method='random',
rp=rp,
trace=TRUE,
rebalance_on=rebalance_period, # uses xts 'endpoints'
trailing_periods=NULL, # calculates from inception
training_period=36) # starts 3 years in to the data history
MRCmETL.RND.w = extractWeights(MRCmETL.RND.t)
MRCmETL.gw = extractGroups(MRCmETL.RND.t)
save(MRCmETL.DE.t,file=paste(resultsdir, 'MRCmETL.DE.t-', Sys.Date(), '-', runname, '.rda',sep=''))
chart.UnStackedBar(MRCmETL.w, rotate="horizontal", colorset=wb13color, space=0, las=2)
MRCmETL=Return.rebalancing(edhec.R, MRCmETL.w)
colnames(MRCmETL) = "MRCmETL"
end_time<-Sys.time()
end_time-start_time
#########################################################################
# Optimization ends here
#########################################################################
braverock/PortfolioAnalytics documentation built on April 18, 2024, 4:09 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
### For Presentation at FactSet's 2013 US Investment Process Symposium
# November 10 - 12 , 2013
# Peter Carl
### Load the necessary packages
# Include optimizer packages
require(PortfolioAnalytics)
require(quantmod)
library(DEoptim)
library(ROI)
require(ROI.plugin.quadprog)
require(ROI.plugin.glpk)
# ... and multi-core packages
require(foreach)
require(doMC)
registerDoMC(6)
# Available on r-forge
# require(FactorAnalytics) # development version > build
### Set script constants
runname='historical.moments'
rebalance_period = 'years' #'quarters' # uses endpoints identifiers from xts; how to do semi-annual?
clean = "boudt" # "none"
permutations = 2000
p=1-1/12 # set confidence for VaR/mETL for monthly data
### Description
# Seven assets, in this case hedge fund indexes representing different styles and two 'types':
# Convertible Arbitrage | Non-directional
# Equity Market Neutral | Non-directional
# Fixed Income Arbitrage | Non-directional
# Event Driven | Non-directional
# CTA Global | Directional
# Global Macro | Directional
# Long/Short Equity | Directional
# see analyze.HFindexes.R for more detail
# Set up seven objectives as buoy portfolios:
# - Equal contribution to...
# 1 Weight
# 2 Variance
# 3 Risk (mETL)
# - Reward to Risk ratio of...
# 4 Mean-Variance
# 5 Mean-mETL
# - Minimum...
# 6 Variance
# 7 Risk (mETL)
# Add constraints
# - Box constraints - 5% to 30%
# - Group constraints - Non-directional constrained to 20-70%; Directional between 10-50%
# - Rebalancing period - quarterly
# - Turnover constraints - TBD
#------------------------------------------------------------------------
# Set up an initial portfolio object with constraints and objectives using
# v2 specification
# Create initial portfolio object used to initialize ALL the bouy portfolios
init.portf <- portfolio.spec(assets=colnames(R),
weight_seq=generatesequence(by=0.005)
)
# Add leverage constraint
init.portf <- add.constraint(portfolio=init.portf,
type="leverage"
)
# Add box constraint
init.portf <- add.constraint(portfolio=init.portf,
type="box",
min=0.05,
max=0.3
)
# Add group constraint
init.portf <- add.constraint(portfolio=init.portf, type="group",
groups=list(c(1:4),
c(5:7)),
group_min=c(0.25,.05),
group_max=c(0.85,0.75)
)
# print(init.portf)
# summary(init.portf)
#------------------------------------------------------------------------
### Construct BUOY 1: Constrained Mean-StdDev Portfolio - using ROI
# Add the return and sd objectives to the constraints created above
MeanSD.portf <- add.objective(portfolio=init.portf,
type="return", # the kind of objective this is
name="mean" # name of the function
)
MeanSD.portf <- add.objective(portfolio=MeanSD.portf,
type="risk", # the kind of objective this is
name="var", # name of the function
arguments=list(clean=clean)
)
### Construct BUOY 2: Constrained Mean-mETL Portfolio - using ROI
#@ Cannot maximize mean return per unit ETL with ROI, consider using
#@ random portfolios or DEoptim. - RB
# Add the return and mETL objectives
MeanmETL.portf <- add.objective(portfolio=init.portf,
type="return", # the kind of objective this is
name="mean" # name of the function
, multiplier=-12
)
MeanmETL.portf <- add.objective(portfolio=MeanmETL.portf,
type="risk", # the kind of objective this is
name="ES", # the function to minimize
arguments=list(p=p, clean=clean)
)
### Construct BUOY 3: Constrained Minimum Variance Portfolio - using ROI
# Add the variance objective
MinSD.portf <- add.objective(portfolio=init.portf,
type="risk", # the kind of objective this is
name="var", # name of the function
arguments=list(p=p, clean=clean)
)
### Construct BUOY 4: Constrained Minimum mETL Portfolio - using ROI
# Add the mETL objective
MinmETL.portf <- add.objective(portfolio=init.portf,
type="risk", # the kind of objective this is
name="ES", # the function to minimize
arguments=list(p=p, clean=clean)
)
### Construct BUOY 5: Constrained Equal Variance Contribution Portfolio - using RP
#@ - Add the sub-objectives first. Adding these 3 objectives means that we are
#@ maximizing mean per unit StdDev with equal volatility contribution portfolios. - RB
# Without a sub-objective, we get a somewhat undefined result, since there are (potentially) many Equal SD contribution portfolios.
# MRCSD.portf <- add.objective(portfolio=init.portf,
# type="risk",
# name="StdDev"
# ) # OR
# MRCSD.portf <- add.objective(portfolio=init.portf,
# type="return",
# name="mean"
# )
MRCSD.portf <- add.objective(portfolio=init.portf,
type="risk_budget",
name="StdDev",
min_concentration=TRUE,
arguments = list(clean=clean)
)
# MRCSD.portf$constraints[[1]]$min_sum = 0.99 # set to speed up RP
# MRCSD.portf$constraints[[1]]$max_sum = 1.01
### Construct BUOY 6: Constrained Equal mETL Contribution Portfolio - using RP
#@ Add the sub-objectives first. These should be added to the MRCmETL portfolio.
#@ All objectives below mean that we are maximizing mean return per unit ES with
#@ equal ES contribution. - RB
# Without a sub-objective, we get a somewhat undefined result, since there are (potentially) many Equal SD contribution portfolios.
# MRCmETL.portf <- add.objective(portfolio=init.portf,
# type="risk",
# name="ES"
# ) # OR
MRCmETL.portf <- add.objective(portfolio=init.portf,
type="return",
name="mean"
)
MRCmETL.portf <- add.objective(MRCmETL.portf,
type="risk_budget",
name="ES",
min_concentration=TRUE,
arguments = list(p=p, clean=clean)
)
# Calculate portfolio variance, but don't use it in the objective; used only for plots
MRCmETL.portf <- add.objective(portfolio=MRCmETL.portf,
type="risk", # the kind of objective this is
name="StdDev", # the function to minimize
enabled=TRUE, # enable or disable the objective
multiplier=0, # calculate it but don't use it in the objective
arguments=list(clean=clean)
)
# MRCmETL.portf$constraints[[1]]$min_sum = 0.99 # set to speed up RP
# MRCmETL.portf$constraints[[1]]$max_sum = 1.01
### Construct BUOY 7: Equal Weight Portfolio
# There's only one, so create a portfolio object with all the objectives we want calculated.
EqWgt.portf <- portfolio.spec(assets=colnames(R))
EqWgt.portf <- add.constraint(portfolio=EqWgt.portf, type="leverage", min_sum=1, max_sum=1)
EqWgt.portf <- add.objective(portfolio=EqWgt.portf, type="return", name="mean")
EqWgt.portf <- add.objective(portfolio=EqWgt.portf, type="risk_budget", name="ES", arguments=list(p=p, clean=clean))
EqWgt.portf <- add.objective(portfolio=EqWgt.portf, type="risk_budget", name="StdDev", arguments=list(clean=clean))
### Construct BUOY 8: Inverse Volatility Portfolio
# There's only one, so create a portfolio object with all the objectives we want calculated.
VolWgt.portf <- portfolio.spec(assets=colnames(R))
VolWgt.portf <- add.constraint(portfolio=VolWgt.portf, type="leverage", min_sum=0.99, max_sum=1.01)
VolWgt.portf <- add.objective(portfolio=VolWgt.portf, type="return", name="mean")
VolWgt.portf <- add.objective(portfolio=VolWgt.portf, type="risk_budget", name="ES", arguments=list(p=p, clean=clean))
VolWgt.portf <- add.objective(portfolio=VolWgt.portf, type="risk_budget", name="StdDev", arguments=list(clean=clean))
# REMOVED - to much to show already
# ### Construct RISK BUDGET Portfolio
# ConstrConcmETL.portf <- portfolio.spec(assets=colnames(R),
# weight_seq=generatesequence(by=0.005)
# )
# # Add leverage constraint
# ConstrConcmETL.portf <- add.constraint(portfolio=RiskBudget.portf,
# type="leverage",
# min_sum=0.99, # set to speed up RP, DE
# max_sum=1.01
# )
# # Establish position bounds
# ConstrConcmETL.portf <- add.constraint(portfolio=ConstrConcmETL.portf,
# type="box",
# min=0.01, # leave relatively unconstrained
# max=1.0
# )
# # Maximize mean return
# ConstrConcmETL.portf <- add.objective(portfolio=ConstrConcmETL.portf,
# type="return", # maximize return
# name="mean",
# multiplier=12
# )
# # Add a risk measure
# # Use ES to be consistent with risk measures in other BUOY portfolios
# ConstrConcmETL.portf <- add.objective(portfolio=ConstrConcmETL.portf,
# type="risk",
# name="ETL", # using a different name to avoid clobbering slot below, workaround for bug
# multiplier=1,
# arguments = list(p=p, clean=clean)
# )
#
# # Set contribution limits
# ConstrConcmETL.portf <- add.objective(portfolio=ConstrConcmETL.portf,
# type="risk_budget",
# name="ES",
# max_prisk=0.3, # Sets the maximum percentage contribution to risk
# arguments = list(p=p, clean=clean)
# )
# Calculate portfolio variance, but don't use it in the objective; used only for plots
# ConstrConcmETL.portf <- add.objective(portfolio=ConstrConcmETL.portf,
# type="risk", # the kind of objective this is
# name="StdDev", # the function to minimize
# enabled=TRUE, # enable or disable the objective
# multiplier=0, # calculate it but don't use it in the objective
# arguments=list(clean=clean)
# )
#------------------------------------------------------------------------
### Evaluate portfolio objective objects
#------------------------------------------------------------------------
# Generate a single set of random portfolios to evaluate against all RP constraint sets
print(paste('constructing random portfolios at',Sys.time()))
# Modify the init.portf specification to get RP running
rp.portf <- init.portf
# rp.portf$constraints[[1]]$min_sum = 0.99 # set to speed up RP
# rp.portf$constraints[[1]]$max_sum = 1.01
rp.portf$constraints[[1]]$min_sum = 1.00 # for more accuracy
rp.portf$constraints[[1]]$max_sum = 1.00
# rp = random_portfolios(portfolio=rp.portf, permutations=30000, max_permutations=400) # will get fewer with less accuracy
load(file=paste(resultsdir,'random-portfolios-2013-10-05.historical.moments.rda'))
rp.mean = apply(rp, 1, function(w) mean(R %*% w))
rp.sd = apply(rp, 1, function(x) StdDev(R=R, weights=x, p=p, clean=clean))
plot(rp.sd, rp.mean, col="darkgray", cex=0.5)
# REMOVED: This was fruitless
# rp1 = random_portfolios(portfolio=rp.portf, permutations=10000, max_permutations=400, rp_method="sample")
# rp1.mean = apply(rp1, 1, function(w) mean(R %*% w))
# rp1.sd = apply(rp1, 1, function(x) StdDev(R=R, weights=x, p=p))
# rp1.etl=NULL; for(i in 1:NROW(rp1)) {rp1.etl[i]=ETL(R=R, weights=as.vector(rp1[i,]), p=p, portfolio_method="component")[[1]]}
# plot(rp1.sd, rp1.mean, col="gray", cex=0.5)
#
# rp2 = random_portfolios(portfolio=rp.portf, permutations=10000, max_permutations=400, rp_method="simplex", fev=2)
# rp2.mean = apply(rp2, 1, function(w) mean(R %*% w))
# rp2.sd = apply(rp2, 1, function(x) StdDev(R=R, weights=x, p=p))
# points(rp2.sd,rp2.mean, col="blue", cex=0.5)
#
# rp3 = random_portfolios(portfolio=rp.portf, permutations=10000, max_permutations=400, rp_method="grid")
# rp3.mean = apply(rp3, 1, function(w) mean(R %*% w))
# rp3.sd = apply(rp3, 1, function(x) StdDev(R=R, weights=x, p=p))
# points(rp3.sd,rp3.mean, col="green", cex=0.5)
# print(paste('done constructing random portfolios at',Sys.time()))
# save(rp,file=paste(resultsdir, 'random-portfolios-', Sys.Date(), '-', runname, '.rda',sep=''))
start_time<-Sys.time()
print(paste('Starting optimization at',Sys.time()))
### Evaluate BUOY 1: Constrained Mean-StdDev Portfolio - with ROI
MeanSD.ROI<-optimize.portfolio(R=R,
portfolio=MeanSD.portf,
optimize_method='ROI',
trace=TRUE
)
plot(MeanSD.ROI, risk.col="StdDev", return.col="mean", rp=permutations, chart.assets=TRUE, main="Mean-Volatility Portfolio")
save(MeanSD.ROI,file=paste(resultsdir, 'MeanSD-', Sys.Date(), '-', runname, '.rda',sep='')) # Save the results
print(paste('Completed meanSD optimization at',Sys.time(),'moving on to meanmETL'))
### Evaluate BUOY 2: Constrained Mean-mETL Portfolio - with DE or RND
# Not possible with ROI - RB
# MeanmETL.ROI<-optimize.portfolio(R=R,
# portfolio=MeanmETL.portf,
# optimize_method='ROI',
# trace=TRUE, verbose=TRUE
# )
#
# So use random portfolios instead
MeanmETL.RND<-optimize.portfolio(R=R,
portfolio=MeanmETL.portf,
optimize_method='random',
rp=rp,
trace=TRUE
)
plot(MeanmETL.RND, risk.col="StdDev", return.col="mean", rp=permutations, chart.assets=TRUE, main="Mean-mETL Portfolio")
plot(MeanmETL.RND, risk.col="ES", return.col="mean", rp=permutations, chart.assets=TRUE, main="Mean-mETL Portfolio")
save(MeanmETL.RND,file=paste(resultsdir, 'MeanETL-RP-', Sys.Date(), '-', runname, '.rda',sep=''))
chart.EfficientFrontier(MeanmETL.RND)
# OR with DE optim
MeanmETL.DE<-optimize.portfolio(R=R,
portfolio=MeanmETL.portf,
optimize_method='DEoptim',
search_size=20000,
initialpop=rp[1:50,] # seed with a starting population that we know fits the constraint space
)
plot(MeanmETL.DE, risk.col="StdDev", return.col="mean", chart.assets=TRUE, main="Mean-mETL Portfolio")
plot(MeanmETL.DE, risk.col="ES", return.col="mean", chart.assets=TRUE, main="Mean-mETL Portfolio")
save(MeanmETL.DE,file=paste(resultsdir, 'MeanETL-DE-', Sys.Date(), '-', runname, '.rda',sep=''))
print(paste('Completed meanmETL optimization at',Sys.time(),'moving on to MinSD'))
### Evaluate BUOY 3: Constrained Minimum Variance Portfolio - with ROI
MinSD.ROI<-optimize.portfolio(R=R,
portfolio=MinSD.portf,
optimize_method='ROI',
trace=TRUE
) #
plot(MinSD.ROI, risk.col="StdDev", return.col="mean", rp=permutations, chart.assets=TRUE, main="Minimum Volatility Portfolio with ROI")
save(MinSD.ROI,file=paste(resultsdir, 'MinSD-', Sys.Date(), '-', runname, '.rda',sep=''))
print(paste('Completed MinSD optimization at',Sys.time(),'moving on to MinmETL'))
# OR with random portfolios
# MinSD.RND<-optimize.portfolio(R=R,
# portfolio=MinSD.portf,
# optimize_method='random',
# rp=rp,
# trace=TRUE, verbose=TRUE
# ) #
# plot(MinSD.RND, risk.col="StdDev", return.col="mean", chart.assets=TRUE, main="Minimum Volatility Portfolio with RP")
### Evaluate BUOY 4: Constrained Minimum mETL Portfolio - with ROI
MinmETL.ROI<-optimize.portfolio(R=R,
portfolio=MinmETL.portf,
optimize_method='ROI',
trace=TRUE, verbose=TRUE
)
plot(MinmETL.ROI, risk.col="StdDev", return.col="mean", rp=permutations, chart.assets=TRUE, main="Minimum mETL Portfolio")
plot(MinmETL.ROI, risk.col="ES", return.col="mean", rp=permutations, chart.assets=TRUE, main="Minimum mETL Portfolio")
save(MinmETL.ROI,file=paste(resultsdir, 'MinmETL-', Sys.Date(), '-', runname, '.rda',sep=''))
print(paste('Completed MinmETL optimization at',Sys.time(),'moving on to MRCSD'))
### Evaluate BUOY 5: Constrained Equal Variance Contribution Portfolio - with RP
# MRCSD.RND<-optimize.portfolio(R=R,
# portfolio=MRCSD.portf,
# optimize_method='random',
# rp=rp,
# trace=TRUE
# )
# plot(MRCSD.RND, risk.col="StdDev", return.col="mean", chart.assets=TRUE, main="Equal Volatility Contribution Portfolio")
# chart.RiskBudget(MRCSD.RND, risk.type="percentage", neighbors=25)
# save(MRCSD.RND,file=paste(resultsdir, 'MRCSD.RND-', Sys.Date(), '-', runname, '.rda',sep=''))
# ... not a very satisfying solution
# OR DE optim - this gets very close (a nice, straight line), so use it
MRCSD.DE<-optimize.portfolio(R=R,
portfolio=MRCSD.portf,
optimize_method='DEoptim',
search_size=20000,
itermax=400,
initialpop=rp[1:50,], # seed with a starting population that we know fits the constraint space
trace=FALSE
)
plot(MRCSD.DE, risk.col="StdDev", return.col="mean", chart.assets=TRUE, main="Equal Volatility Contribution Portfolio")
chart.RiskBudget(MRCSD.DE, risk.type="percentage", neighbors=25)
save(MRCSD.DE,file=paste(resultsdir, 'MRCSD.DE-', Sys.Date(), '-', runname, '.rda',sep=''))
print(paste('Completed MRCSD optimization at',Sys.time(),'moving on to MRCmETL'))
### Evaluate BUOY 6: Constrained Equal mETL Contribution Portfolio - with RP
MRCmETL.RND<-optimize.portfolio(R=R,
portfolio=MRCmETL.portf,
optimize_method='random',
rp=rp,
trace=TRUE
) #
plot(MRCmETL.RND, risk.col="StdDev", return.col="mean", chart.assets=TRUE, main="Equal mETL Contribution Portfolio")
plot(MRCmETL.RND, risk.col="ES", return.col="mean", chart.assets=TRUE, main="Equal mETL Contribution Portfolio")
chart.RiskBudget(MRCmETL.RND, neighbors=25)
save(MRCmETL.RND,file=paste(resultsdir, 'MRCmETL-', Sys.Date(), '-', runname, '.rda',sep=''))
# OR DE optim -
MRCmETL.DE<-optimize.portfolio(R=R,
portfolio=MRCmETL.portf,
optimize_method='DEoptim',
search_size=20000,
NP=200,
initialpop=rp[1:50,], # seed with a starting population that we know fits the constraint space
trace=FALSE
)
plot(MRCmETL.DE, risk.col="StdDev", return.col="mean", chart.assets=TRUE, main="Equal Volatility Contribution Portfolio")
chart.RiskBudget(MRCmETL.DE, risk.type="percentage", neighbors=25)
save(MRCmETL.DE,file=paste(resultsdir, 'MRCmETL.DE-', Sys.Date(), '-', runname, '.rda',sep=''))
# # test it unconstrained:
# unconstr.portf <- portfolio.spec(assets=colnames(R),
# weight_seq=generatesequence(by=0.005)
# )
# unconstr.portf <- add.constraint(portfolio=unconstr.portf,
# type="leverage",
# min_sum=0.99, # set to speed up RP
# max_sum=1.01
# )
# # Establish position bounds
# unconstr.portf <- add.constraint(portfolio=unconstr.portf,
# type="box",
# min=0.01,
# max=1.0
# )
# MRCmETLun.portf <- add.objective(portfolio=unconstr.portf,
# type="return",
# name="mean"
# )
# MRCmETLun.portf <- add.objective(MRCmETL.portf,
# type="risk_budget",
# name="ES",
# min_concentration=TRUE,
# arguments = list(p=p, clean=clean)
# )
#
# # ...in DE optim -
# MRCmETLun.DE<-optimize.portfolio(R=R,
# portfolio=MRCmETLun.portf,
# optimize_method='DEoptim',
# search_size=20000,
# NP=200,
# initialpop=rp[1:50,], # seed with a starting population that we know fits the constraint space
# trace=FALSE
# )
print(paste('Completed MRCmETL optimization at',Sys.time(),'moving on to RiskBudget'))
### Evaluate BUOY 7: Equal Weight Portfolio
# Calculate the objective measures for the equal weight portfolio
EqWgt.opt <- equal.weight(R=R, portfolio=EqWgt.portf)
### Evaluate BUOY 8: Inverse Volatility Portfolio
volatility.weight <- function (R, portfolio, ...)
{
if (!is.portfolio(portfolio))
stop("portfolio object passed in must be of class 'portfolio'")
assets <- portfolio$assets
nassets <- length(assets)
weights <- as.vector((1/StdDev(R))/sum(1/StdDev(R)))
names(weights) <- names(assets)
if (ncol(R) != nassets) {
if (ncol(R) > nassets) {
R <- R[, 1:nassets]
warning("number of assets is less than number of columns in returns object, subsetting returns object.")
}
else {
stop("number of assets is greater than number of columns in returns object")
}
}
out <- constrained_objective(w = weights, R = R, portfolio = portfolio,
trace = TRUE, ...)$objective_measures
return(structure(list(R = R, weights = weights, objective_measures = out,
call = match.call(), portfolio = portfolio), class = c("optimize.portfolio.invol",
"optimize.portfolio")))
}
# Calculate the objective measures for the vol weight portfolio
VolWgt.opt <- volatility.weight(R=R, portfolio=VolWgt.portf)
# REMOVED
# ### Evaluate Constrained Concentration to mETL Portfolio - with DE
# # registerDoSEQ() # turn off parallelization to keep the trace data
# ConstrConcmETL.DE<-optimize.portfolio(R=R,
# portfolio=ConstrConcmETL.portf,
# optimize_method='DEoptim',
# search_size=40000,
# NP=4000,
# itermax=400,
# trace=FALSE
# )
# # list(c=0.25, # speed of crossover adaption (0,1]
# # CR=0.75) # crossover probability [0,1]
# plot(ConstrConcmETL.DE, risk.col="StdDev", return.col="mean")
# plot(ConstrConcmETL.DE, risk.col="ES", return.col="mean") # several outlier portfolios
# chart.RiskBudget(ConstrConcmETL.DE)
# chart.RiskBudget(ConstrConcmETL.DE, risk.type="percentage")
#
# save(ConstrConcmETL.DE,file=paste(resultsdir, 'ConstrConcmETL-', Sys.Date(), '-', runname, '.rda',sep=''))
# print(ConstrConcmETL.DE$elapsed_time)
print('Done with optimizations.')
#------------------------------------------------------------------------
### Extract data from optimizations for analysis
#------------------------------------------------------------------------
# Combine optimization objects
buoys <- combine.optimizations(list(MeanSD=MeanSD.ROI, MeanmETL=MeanmETL.RND, MinSD=MinSD.ROI, MinmETL=MinmETL.ROI, MRCSD=MRCSD.DE, MRCmETL=MRCmETL.DE, VolWgt=VolWgt.opt, EqWgt=EqWgt.opt))
# chart.Weights(buoys, plot.type="bar", ylim=c(0,1))
#
# #@ Chart the portfolios that have mean and ES as objective measures. - RB
# chart.RiskReward(buoys, risk.col="ES")
# #@ Chart the portfolios that have mean and StdDev as objective measures. - RB
# chart.RiskReward(buoys, risk.col="StdDev")
#
# #@ The MRCmETL and RB optimizations would be good to compare because they are
# #@ similar in that they both include component ES as an objective. - RB
# buoyETL <- combine.optimizations(list(MRCmETL=MRCmETL.RND, RB=RiskBudget.DE, EqWgt=EqWgt.opt))
# chart.RiskBudget(buoyETL, match.col="ES", risk.type="percentage", legend.loc="topright")
#
# #@ Compare the equal weight portfolio and the equal SD contribution portfolio. - RB
# buoyStdDev <- combine.optimizations(list(MRCSD=MRCSD.RND, EqWgt=EqWgt.opt))
# chart.RiskBudget(buoyStdDev, match.col="StdDev", risk.type="absolute", legend.loc="topleft")
Wgts = extractWeights(buoys)
### Extract portfolio measures from each objective
# We can't just extract them, because they aren't all calculated
# so fill them in...
buoys.portfmeas = buoys.contrib.sd = buoys.contrib.es = buoys.perc.sd = buoys.perc.es = NULL
for(i in 1:NROW(Wgts)){
mean = sum(colMeans(R)*Wgts[i,])
sd = StdDev(R, weights=Wgts[i,], portfolio_method="component", clean=clean)
es = ES(R, weights=Wgts[i,], method="modified", portfolio_method="component", p=p, clean=clean)
buoys.portfmeas=rbind(buoys.portfmeas, c(mean, sd[[1]][1], es[[1]][1]))
buoys.contrib.sd= rbind(buoys.contrib.sd,sd[[2]])
buoys.contrib.es= rbind(buoys.contrib.es,es[[2]])
buoys.perc.sd = rbind(buoys.perc.sd,sd[[3]])
buoys.perc.es = rbind(buoys.perc.es,es[[3]])
}
colnames(buoys.portfmeas)=c("Mean", "StdDev", "mETL")
rownames(buoys.portfmeas)=
rownames(buoys.contrib.sd)=
rownames(buoys.contrib.es)=
rownames(buoys.perc.sd) =
rownames(buoys.perc.es) = rownames(Wgts)
colnames(buoys.contrib.sd)=
colnames(buoys.contrib.es)=
colnames(buoys.perc.sd) =
colnames(buoys.perc.es) = colnames(Wgts)
# get the RP portfolios with risk and return pre-calculated
xtract = extractStats(MRCmETL.RND)
save(xtract,file=paste(resultsdir, 'xtract-RPs-', Sys.Date(), '-', runname, '.rda',sep=''))
# columnnames = colnames(xtract)
results.names=rownames(buoys.portfmeas)
# by Asset metrics
assets.portfmeas=as.matrix(scatterFUN(R, FUN="mean"))
assets.portfmeas=cbind(assets.portfmeas, scatterFUN(R, FUN="StdDev", clean=clean))
assets.portfmeas=cbind(assets.portfmeas, scatterFUN(R, FUN="ES", clean=clean))
colnames(assets.portfmeas)=c("Mean", "StdDev", "mETL")
rownames(assets.portfmeas)=colnames(Wgts)
#------------------------------------------------------------------------
# Run select buoy optimizations through time
#------------------------------------------------------------------------
#
# Equal Weight
dates=index(R[endpoints(R, on=rebalance_period)])
EqWgt.w = xts(matrix(rep(1/NCOL(R),length(dates)*NCOL(R)), ncol=NCOL(R)), order.by=dates)
colnames(EqWgt.w)= colnames(R)
EqWgt.R=Return.rebalancing(R, EqWgt.w)
chart.StackedBar(EqWgt.w, colorset=wb13color, gap=0)
VolWgt.w = NULL
for(i in 3:length(dates)){
x = volatility.weight(R=R[paste0("::",dates[i]),], portfolio=VolWgt.portf)
VolWgt.w = rbind(VolWgt.w, x$weights)
}
VolWgt.w = as.xts(VolWgt.w, order.by=dates[-1:-2])
VolWgt.R=Return.rebalancing(R, VolWgt.w)
# Equal SD
MRCSD.DE.t = optimize.portfolio.rebalancing(R=R,
portfolio=MRCSD.portf,
optimize_method='DEoptim',
search_size=2000,
NP=200,
initialpop=rp[1:50,], # seed with a starting population that we know fits the constraint space
trace=FALSE,
rebalance_on=rebalance_period, # uses xts 'endpoints'
trailing_periods=NULL, # calculates from inception
training_period=36) # starts 3 years in to the data history
MRCSD.w = extractWeights(MRCSD.DE.t)
MRCSD.gw = extractGroups(MRCSD.DE.t)
save(MRCSD.DE.t,file=paste(resultsdir, 'MRCSD.DE.t-', Sys.Date(), '-', runname, '.rda',sep=''))
chart.UnStackedBar(MRCSD.w, rotate="horizontal", colorset=wb13color, space=0, las=2)
MRCSD.R=Return.rebalancing(edhec.R, MRCSD.w)
colnames(MRCSD) = "MRCSD"
# Extract perc contrib of mES from results object
x=NULL
for(i in 1:length(names(MRCSD.DE.t))) {
x = rbind(x,MRCSD.DE.t[[i]][["objective_measures"]]$StdDev$pct_contrib_StdDev)
}
MRCSD.DE.pct_contrib_StdDev.t = as.xts(x, order.by=as.POSIXct(names(MRCSD.DE.t)))
chart.UnStackedBar(x.xts, rotate="horizontal", colorset=wb13color, space=0, las=2)
# MRC mETL
MRCmETL.DE.t = optimize.portfolio.rebalancing(R=R,
portfolio=MRCmETL.portf,
optimize_method='DEoptim',
search_size=20000,
NP=200,
initialpop=rp[1:50,], # seed with a starting population that we know fits the constraint space
trace=FALSE,
rebalance_on=rebalance_period, # uses xts 'endpoints'
trailing_periods=NULL, # calculates from inception
training_period=36) # starts 3 years in to the data history
MRCmETL.w = extractWeights(MRCmETL.DE.t)
MRCmETL.gw = extractGroups(MRCmETL.DE.t)
save(MRCmETL.DE.t,file=paste(resultsdir, 'MRCmETL.DE.t-', Sys.Date(), '-', runname, '.rda',sep=''))
chart.UnStackedBar(MRCmETL.w, rotate="horizontal", colorset=wb13color, space=0, las=2)
MRCmETL=Return.rebalancing(edhec.R, MRCmETL.w)
colnames(MRCmETL) = "MRCmETL"
MRCmETL.RND.t = optimize.portfolio.rebalancing(R=R,
portfolio=MRCmETL.portf,
optimize_method='random',
rp=rp,
trace=TRUE,
rebalance_on=rebalance_period, # uses xts 'endpoints'
trailing_periods=NULL, # calculates from inception
training_period=36) # starts 3 years in to the data history
MRCmETL.RND.w = extractWeights(MRCmETL.RND.t)
MRCmETL.gw = extractGroups(MRCmETL.RND.t)
save(MRCmETL.DE.t,file=paste(resultsdir, 'MRCmETL.DE.t-', Sys.Date(), '-', runname, '.rda',sep=''))
chart.UnStackedBar(MRCmETL.w, rotate="horizontal", colorset=wb13color, space=0, las=2)
MRCmETL=Return.rebalancing(edhec.R, MRCmETL.w)
colnames(MRCmETL) = "MRCmETL"
end_time<-Sys.time()
end_time-start_time
#########################################################################
# Optimization ends here
#########################################################################
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.