###############################################################################
# R (http://r-project.org/) Numeric Methods for Optimization of Portfolios
#
# Copyright (c) 2004-2014 Brian G. Peterson, Peter Carl, Ross Bennett, Kris Boudt
#
# This library is distributed under the terms of the GNU Public License (GPL)
# for full details see the file COPYING
#
# $Id$
#
###############################################################################
#' @rdname optimize.portfolio
#' @name optimize.portfolio
#' @export
optimize.portfolio_v1 <- function(
R,
constraints,
optimize_method=c("DEoptim","random","ROI","ROI_old","pso","GenSA"),
search_size=20000,
trace=FALSE, ...,
rp=NULL,
momentFUN='set.portfolio.moments_v1'
)
{
optimize_method=optimize_method[1]
tmptrace=NULL
start_t<-Sys.time()
#store the call for later
call <- match.call()
if (is.null(constraints) | !is.constraint(constraints)){
stop("you must pass in an object of class constraints to control the optimization")
}
R <- checkData(R)
N = length(constraints$assets)
if (ncol(R)>N) {
R=R[,names(constraints$assets)]
}
T = nrow(R)
out=list()
weights=NULL
dotargs <-list(...)
# set portfolio moments only once
if(!is.function(momentFUN)){
momentFUN<-match.fun(momentFUN)
}
# TODO FIXME should match formals later
#dotargs <- set.portfolio.moments(R, constraints, momentargs=dotargs)
.mformals <- dotargs
.mformals$R <- R
.mformals$constraints <- constraints
mout <- try((do.call(momentFUN,.mformals)) ,silent=TRUE)
if(inherits(mout,"try-error")) {
message(paste("portfolio moment function failed with message",mout))
} else {
dotargs <- c(dotargs,mout)
}
normalize_weights <- function(weights){
# normalize results if necessary
if(!is.null(constraints$min_sum) | !is.null(constraints$max_sum)){
# the user has passed in either min_sum or max_sum constraints for the portfolio, or both.
# we'll normalize the weights passed in to whichever boundary condition has been violated
# NOTE: this means that the weights produced by a numeric optimization algorithm like DEoptim
# might violate your constraints, so you'd need to renormalize them after optimizing
# we'll create functions for that so the user is less likely to mess it up.
##' NOTE: need to normalize in the optimization wrapper too before we return, since we've normalized in here
##' In Kris' original function, this was manifested as a full investment constraint
if(!is.null(constraints$max_sum) & constraints$max_sum != Inf ) {
max_sum=constraints$max_sum
if(sum(weights)>max_sum) { weights<-(max_sum/sum(weights))*weights } # normalize to max_sum
}
if(!is.null(constraints$min_sum) & constraints$min_sum != -Inf ) {
min_sum=constraints$min_sum
if(sum(weights)<min_sum) { weights<-(min_sum/sum(weights))*weights } # normalize to min_sum
}
} # end min_sum and max_sum normalization
return(weights)
}
if(optimize_method=="DEoptim"){
stopifnot("package:DEoptim" %in% search() || require("DEoptim",quietly = TRUE) )
# DEoptim does 200 generations by default, so lets set the size of each generation to search_size/200)
if(hasArg(itermax)) itermax=match.call(expand.dots=TRUE)$itermax else itermax=N*50
NP = round(search_size/itermax)
if(NP<(N*10)) NP <- N*10
if(NP>2000) NP=2000
if(!hasArg(itermax)) {
itermax<-round(search_size/NP)
if(itermax<50) itermax=50 #set minimum number of generations
}
#check to see whether we need to disable foreach for parallel optimization, esp if called from inside foreach
if(hasArg(parallel)) parallel=match.call(expand.dots=TRUE)$parallel else parallel=TRUE
if(!isTRUE(parallel) && 'package:foreach' %in% search()){
registerDoSEQ()
}
DEcformals <- formals(DEoptim.control)
DEcargs <- names(DEcformals)
if( is.list(dotargs) ){
pm <- pmatch(names(dotargs), DEcargs, nomatch = 0L)
names(dotargs[pm > 0L]) <- DEcargs[pm]
DEcformals$NP <- NP
DEcformals$itermax <- itermax
DEcformals[pm] <- dotargs[pm > 0L]
if(!hasArg(strategy)) DEcformals$strategy=6 # use DE/current-to-p-best/1
if(!hasArg(reltol)) DEcformals$reltol=.000001 # 1/1000 of 1% change in objective is significant
if(!hasArg(steptol)) DEcformals$steptol=round(N*1.5) # number of assets times 1.5 tries to improve
if(!hasArg(c)) DEcformals$c=.4 # JADE mutation parameter, this could maybe use some adjustment
if(!hasArg(storepopfrom)) DEcformals$storepopfrom=1
if(isTRUE(parallel) && 'package:foreach' %in% search()){
if(!hasArg(parallelType) ) DEcformals$parallelType='auto' #use all cores
if(!hasArg(packages) ) DEcformals$packages <- names(sessionInfo()$otherPkgs) #use all packages
}
#TODO FIXME also check for a passed in controlDE list, including checking its class, and match formals
}
if(isTRUE(trace)) {
#we can't pass trace=TRUE into constrained objective with DEoptim, because it expects a single numeric return
tmptrace=trace
assign('.objectivestorage', list(), pos='.GlobalEnv')
trace=FALSE
}
# get upper and lower weights parameters from constraints
upper = constraints$max
lower = constraints$min
if(hasArg(rpseed)){
seed <- match.call(expand.dots=TRUE)$rpseed
DEcformals$initialpop <- seed
rpseed <- FALSE
} else {
rpseed <- TRUE
}
if(hasArg(rpseed) & isTRUE(rpseed)) {
# initial seed population is generated with random_portfolios function
if(hasArg(eps)) eps=match.call(expand.dots=TRUE)$eps else eps = 0.01
rpconstraint<-constraint(assets=length(lower), min_sum=constraints$min_sum-eps, max_sum=constraints$max_sum+eps,
min=lower, max=upper, weight_seq=generatesequence())
rp <- random_portfolios_v1(rpconstraints=rpconstraint,permutations=NP)
DEcformals$initialpop=rp
}
controlDE <- do.call(DEoptim.control,DEcformals)
# minw = try(DEoptim( constrained_objective , lower = lower[1:N] , upper = upper[1:N] , control = controlDE, R=R, constraints=constraints, ...=...)) # add ,silent=TRUE here?
minw = try(DEoptim( constrained_objective_v1 , lower = lower[1:N] , upper = upper[1:N] , control = controlDE, R=R, constraints=constraints, nargs = dotargs , ...=...)) # add ,silent=TRUE here?
if(inherits(minw,"try-error")) { minw=NULL }
if(is.null(minw)){
message(paste("Optimizer was unable to find a solution for target"))
return (paste("Optimizer was unable to find a solution for target"))
}
if(isTRUE(tmptrace)) trace <- tmptrace
weights = as.vector( minw$optim$bestmem)
weights <- normalize_weights(weights)
names(weights) = colnames(R)
out = list(weights=weights, objective_measures=constrained_objective_v1(w=weights,R=R,constraints,trace=TRUE)$objective_measures,out=minw$optim$bestval, call=call)
if (isTRUE(trace)){
out$DEoutput=minw
out$DEoptim_objective_results<-try(get('.objectivestorage',pos='.GlobalEnv'),silent=TRUE)
rm('.objectivestorage',pos='.GlobalEnv')
}
} ## end case for DEoptim
if(optimize_method=="random"){
#' call random_portfolios() with constraints and search_size to create matrix of portfolios
if(missing(rp) | is.null(rp)){
rp<-random_portfolios_v1(rpconstraints=constraints,permutations=search_size)
}
#' store matrix in out if trace=TRUE
if (isTRUE(trace)) out$random_portfolios<-rp
#' write foreach loop to call constrained_objective() with each portfolio
if ("package:foreach" %in% search() & !hasArg(parallel)){
rp_objective_results<-foreach(ii=1:nrow(rp), .errorhandling='pass') %dopar% constrained_objective_v1(w=rp[ii,],R,constraints,trace=trace,...=dotargs)
} else {
rp_objective_results<-apply(rp, 1, constrained_objective_v1, R=R, constraints=constraints, trace=trace, ...=dotargs)
}
#' if trace=TRUE , store results of foreach in out$random_results
if(isTRUE(trace)) out$random_portfolio_objective_results<-rp_objective_results
#' loop through results keeping track of the minimum value of rp_objective_results[[objective]]$out
search<-vector(length=length(rp_objective_results))
# first we construct the vector of results
for (i in 1:length(search)) {
if (isTRUE(trace)) {
search[i]<-ifelse(try(rp_objective_results[[i]]$out),rp_objective_results[[i]]$out,1e6)
} else {
search[i]<-as.numeric(rp_objective_results[[i]])
}
}
# now find the weights that correspond to the minimum score from the constrained objective
# and normalize_weights so that we meet our min_sum/max_sum constraints
if (isTRUE(trace)) {
min_objective_weights<- try(normalize_weights(rp_objective_results[[which.min(search)]]$weights))
} else {
min_objective_weights<- try(normalize_weights(rp[which.min(search),]))
}
#' re-call constrained_objective on the best portfolio, as above in DEoptim, with trace=TRUE to get results for out list
out$weights<-min_objective_weights
out$objective_measures<-try(constrained_objective_v1(w=min_objective_weights,R=R,constraints,trace=TRUE)$objective_measures)
out$call<-call
#' construct out list to be as similar as possible to DEoptim list, within reason
} ## end case for random
if(optimize_method == "ROI_old"){
# This will take a new constraint object that is of the same structure of a
# ROI constraint object, but with an additional solver arg.
# then we can do something like this
print("ROI_old is going to be depricated.")
roi.result <- ROI::ROI_solve(x=constraints$constrainted_objective, constraints$solver)
weights <- roi.result$solution
names(weights) <- colnames(R)
out$weights <- weights
out$objective_measures <- roi.result$objval
out$call <- call
} ## end case for ROI_old
if(optimize_method == "ROI"){
# This takes in a regular constraint object and extracts the desired business objectives
# and converts them to matrix form to be inputed into a closed form solver
# Applying box constraints
bnds <- list(lower = list(ind = seq.int(1L, N), val = as.numeric(constraints$min)),
upper = list(ind = seq.int(1L, N), val = as.numeric(constraints$max)))
# retrive the objectives to minimize, these should either be "var" and/or "mean"
# we can eight miniminze variance or maximize quiadratic utility (we will be minimizing the neg. quad. utility)
moments <- list(mean=rep(0, N))
alpha <- 0.05
target <- NA
for(objective in constraints$objectives){
if(objective$enabled){
if(!any(c(objective$name == "mean", objective$name == "var", objective$name == "CVaR")))
stop("ROI only solves mean, var, or sample CVaR type business objectives, choose a different optimize_method.")
# I'm not sure what changed, but moments$mean used to be a vector of the column means
# now it is a scalar value of the mean of the entire R object
if(objective$name == "mean"){
moments[[objective$name]] <- try(as.vector(apply(R, 2, "mean", na.rm=TRUE)), silent=TRUE)
} else {
moments[[objective$name]] <- try(eval(as.symbol(objective$name))(R), silent=TRUE)
}
target <- ifelse(!is.null(objective$target),objective$target, target)
alpha <- ifelse(!is.null(objective$alpha), objective$alpha, alpha)
lambda <- ifelse(!is.null(objective$risk_aversion), objective$risk_aversion, 1)
}
}
plugin <- ifelse(any(names(moments)=="var"), "quadprog", "glpk")
if(plugin == "quadprog") ROI_objective <- ROI::Q_objective(Q=2*lambda*moments$var, L=-moments$mean)
if(plugin == "glpk") ROI_objective <- ROI::L_objective(L=-moments$mean)
Amat <- rbind(rep(1, N), rep(1, N))
dir.vec <- c(">=","<=")
rhs.vec <- c(constraints$min_sum, constraints$max_sum)
if(!is.na(target)) {
Amat <- rbind(Amat, moments$mean)
dir.vec <- c(dir.vec, "==")
rhs.vec <- c(rhs.vec, target)
}
if(try(!is.null(constraints$groups), silent=TRUE)){
if(sum(constraints$groups) != N)
stop("Number of assets in each group needs to sum to number of total assets.")
n.groups <- length(constraints$groups)
if(!all(c(length(constraints$cLO),length(constraints$cLO)) == n.groups) )
stop("Number of group constraints exceeds number of groups.")
Amat.group <- matrix(0, nrow=n.groups, ncol=N)
k <- 1
l <- 0
for(i in 1:n.groups){
j <- constraints$groups[i]
Amat.group[i, k:(l+j)] <- 1
k <- l + j + 1
l <- k - 1
}
if(is.null(constraints$cLO)) cLO <- rep(-Inf, n.groups)
if(is.null(constraints$cUP)) cUP <- rep(Inf, n.groups)
Amat <- rbind(Amat, Amat.group, -Amat.group)
dir.vec <- c(dir.vec, rep(">=", (n.groups + n.groups)))
rhs.vec <- c(rhs.vec, constraints$cLO, -constraints$cUP)
}
if(any(names(moments)=="CVaR")) {
Rmin <- ifelse(is.na(target), 0, target)
ROI_objective <- ROI::L_objective(c(rep(0,N), rep(1/(alpha*T),T), 1))
Amat <- cbind(rbind(1, 1, moments$mean, coredata(R)), rbind(0, 0, 0, cbind(diag(T), 1)))
dir.vec <- c(">=","<=",">=",rep(">=",T))
rhs.vec <- c(constraints$min_sum, constraints$max_sum, Rmin ,rep(0, T))
if(try(!is.null(constraints$groups), silent=TRUE)){
zeros <- matrix(0, nrow=n.groups, ncol=(T+1))
Amat <- rbind(Amat, cbind(Amat.group, zeros), cbind(-Amat.group, zeros))
dir.vec <- c(dir.vec, rep(">=", (n.groups + n.groups)))
rhs.vec <- c(rhs.vec, constraints$cLO, -constraints$cUP)
}
}
opt.prob <- ROI::OP(objective=ROI_objective,
constraints=ROI::L_constraint(L=Amat, dir=dir.vec, rhs=rhs.vec),
bounds=bnds)
roi.result <- ROI::ROI_solve(x=opt.prob, solver=plugin)
weights <- roi.result$solution[1:N]
names(weights) <- colnames(R)
out$weights <- weights
out$out <- roi.result$objval
out$call <- call
} ## end case for ROI
## case if method=pso---particle swarm
if(optimize_method=="pso"){
stopifnot("package:pso" %in% search() || require("pso",quietly = TRUE) )
if(hasArg(maxit)) maxit=match.call(expand.dots=TRUE)$maxit else maxit=N*50
controlPSO <- list(trace=FALSE, fnscale=1, maxit=1000, maxf=Inf, abstol=-Inf, reltol=0)
PSOcargs <- names(controlPSO)
if( is.list(dotargs) ){
pm <- pmatch(names(dotargs), PSOcargs, nomatch = 0L)
names(dotargs[pm > 0L]) <- PSOcargs[pm]
controlPSO$maxit <- maxit
controlPSO[pm] <- dotargs[pm > 0L]
if(!hasArg(reltol)) controlPSO$reltol <- .0001 # 1/100 of 1% change in objective is insignificant enough to restart a swarm
#NOTE reltol has a different meaning for pso than it has for DEoptim. for DEoptim, reltol is a stopping criteria, for pso,
# it is a restart criteria.
if(!hasArg(s)) controlPSO$s<-N*10 #swarm size
if(!hasArg(maxit.stagnate)) controlPSO$maxit.stagnate <- controlPSO$s #stopping criteria
if(hasArg(trace) && try(trace==TRUE,silent=TRUE)) controlPSO$trace <- TRUE
if(hasArg(trace) && isTRUE(trace)) {
controlPSO$trace <- TRUE
controlPSO$trace.stats=TRUE
}
}
# get upper and lower weights parameters from constraints
upper <- constraints$max
lower <- constraints$min
minw = try(psoptim( par = rep(NA, N), fn = constrained_objective_v1 , R=R, constraints=constraints,
lower = lower[1:N] , upper = upper[1:N] , control = controlPSO)) # add ,silent=TRUE here?
if(inherits(minw,"try-error")) { minw=NULL }
if(is.null(minw)){
message(paste("Optimizer was unable to find a solution for target"))
return (paste("Optimizer was unable to find a solution for target"))
}
weights <- as.vector( minw$par)
weights <- normalize_weights(weights)
names(weights) <- colnames(R)
out = list(weights=weights,
objective_measures=constrained_objective_v1(w=weights,R=R,constraints,trace=TRUE)$objective_measures,
out=minw$value,
call=call)
if (isTRUE(trace)){
out$PSOoutput=minw
}
} ## end case for pso
## case if method=GenSA---Generalized Simulated Annealing
if(optimize_method=="GenSA"){
stopifnot("package:GenSA" %in% search() || require("GenSA",quietly = TRUE) )
if(hasArg(maxit)) maxit=match.call(expand.dots=TRUE)$maxit else maxit=N*50
controlGenSA <- list(maxit = 5000, threshold.stop = NULL, temp = 5230,
visiting.param = 2.62, acceptance.param = -5, max.time = NULL,
nb.stop.improvement = 1e+06, smooth = TRUE, max.call = 1e+07,
verbose = FALSE)
GenSAcargs <- names(controlGenSA)
if( is.list(dotargs) ){
pm <- pmatch(names(dotargs), GenSAcargs, nomatch = 0L)
names(dotargs[pm > 0L]) <- GenSAcargs[pm]
controlGenSA$maxit <- maxit
controlGenSA[pm] <- dotargs[pm > 0L]
if(hasArg(trace) && try(trace==TRUE,silent=TRUE)) controlGenSA$verbose <- TRUE
}
upper <- constraints$max
lower <- constraints$min
minw = try(GenSA( par = rep(1/N, N), lower = lower[1:N] , upper = upper[1:N], control = controlGenSA,
fn = constrained_objective_v1 , R=R, constraints=constraints)) # add ,silent=TRUE here?
if(inherits(minw,"try-error")) { minw=NULL }
if(is.null(minw)){
message(paste("Optimizer was unable to find a solution for target"))
return (paste("Optimizer was unable to find a solution for target"))
}
weights <- as.vector(minw$par)
weights <- normalize_weights(weights)
names(weights) <- colnames(R)
out = list(weights=weights,
objective_measures=constrained_objective_v1(w=weights,R=R,constraints,trace=TRUE)$objective_measures,
out=minw$value,
call=call)
if (isTRUE(trace)){
out$GenSAoutput=minw
}
} ## end case for GenSA
end_t<-Sys.time()
# print(c("elapsed time:",round(end_t-start_t,2),":diff:",round(diff,2), ":stats: ", round(out$stats,4), ":targets:",out$targets))
message(c("elapsed time:",end_t-start_t))
out$constraints<-constraints
out$data_summary<-list(first=first(R),last=last(R))
out$elapsed_time<-end_t-start_t
out$end_t<-as.character(Sys.time())
class(out)<-c(paste("optimize.portfolio",optimize_method,sep='.'),"optimize.portfolio")
return(out)
}
##### version 2 of optimize.portfolio #####
#' @rdname optimize.portfolio
#' @export
optimize.portfolio_v2 <- function(
R,
portfolio=NULL,
constraints=NULL,
objectives=NULL,
optimize_method=c("DEoptim","random","ROI","pso","GenSA"),
search_size=20000,
trace=FALSE, ...,
rp=NULL,
momentFUN='set.portfolio.moments',
message=FALSE
)
{
# This is the case where the user has passed in a list of portfolio objects
# for the portfolio argument.
# Loop through the portfolio list and recursively call optimize.portfolio
# Note that I return at the end of this block. I know it is not good practice
# to return before the end of a function, but I am not sure of another way
# to handle a list of portfolio objects with the recursive call to
# optimize.portfolio.
if(inherits(portfolio, "portfolio.list")){
n.portf <- length(portfolio)
opt.list <- vector("list", n.portf)
for(i in 1:length(opt.list)){
if(message) cat("Starting optimization of portfolio ", i, "\n")
opt.list[[i]] <- optimize.portfolio(R=R,
portfolio=portfolio[[i]],
constraints=constraints,
objectives=objectives,
optimize_method=optimize_method,
search_size=search_size,
trace=trace,
...=...,
rp=rp,
momentFUN=momentFUN,
message=message)
}
out <- combine.optimizations(opt.list)
##### return here for portfolio.list because this is a recursive call
##### for optimize.portfolio
return(out)
}
# Detect regime switching portfolio
if(inherits(portfolio, "regime.portfolios")){
regime.switching <- TRUE
regime <- portfolio$regime
if(index(last(R)) %in% index(regime)){
regime.idx <- as.numeric(regime[index(last(R))])[1]
portfolio <- portfolio$portfolio.list[[regime.idx]]
#cat("regime: ", regime.idx, "\n")
} else {
warning("Dates in regime and R do not match, defaulting to first portfolio")
regime.idx <- 1
portfolio <- portfolio$portfolio.list[[regime.idx]]
}
} else {
regime.switching <- FALSE
}
optimize_method <- optimize_method[1]
tmptrace <- NULL
start_t <- Sys.time()
#store the call for later
call <- match.call()
if (!is.null(portfolio) & !is.portfolio(portfolio)){
stop("you must pass in an object of class 'portfolio' to control the optimization")
}
# Check for constraints and objectives passed in separately outside of the portfolio object
if(!is.null(constraints)){
if(inherits(constraints, "v1_constraint")){
if(is.null(portfolio)){
# If the user has not passed in a portfolio, we will create one for them
tmp_portf <- portfolio.spec(assets=constraints$assets)
}
message("constraint object passed in is a 'v1_constraint' object, updating to v2 specification")
portfolio <- update_constraint_v1tov2(portfolio=tmp_portf, v1_constraint=constraints)
# print.default(portfolio)
}
if(!inherits(constraints, "v1_constraint")){
# Insert the constraints into the portfolio object
portfolio <- insert_constraints(portfolio=portfolio, constraints=constraints)
}
}
if(!is.null(objectives)){
# Insert the objectives into the portfolio object
portfolio <- insert_objectives(portfolio=portfolio, objectives=objectives)
}
R <- checkData(R)
N <- length(portfolio$assets)
if (ncol(R) > N) {
R <- R[,names(portfolio$assets)]
}
T <- nrow(R)
# Initialize an empty list used as the return object
out <- list()
weights <- NULL
# Get the constraints from the portfolio object
constraints <- get_constraints(portfolio)
# set portfolio moments only once
# For set.portfolio.moments, we are passing the returns,
# portfolio object, and dotargs. dotargs is a list of arguments
# that are passed in as dots in optimize.portfolio. This was
# causing errors if clean="boudt" was specified in an objective
# and an argument such as itermax was passed in as dots to
# optimize.portfolio. See r2931
if(!is.function(momentFUN)){
momentFUN <- match.fun(momentFUN)
}
# **
# When an ES/ETL/CVaR problem is being solved by a linear solver, the higher
# moments do not need to be calculated. The moments are very compute
# intensive and slow down the optimization problem.
# match the args for momentFUN
.formals <- formals(momentFUN)
.formals <- modify.args(formals=.formals, arglist=NULL, ..., dots=FALSE)
# ** pass ROI=TRUE to set.portfolio.moments so the moments are not calculated
if(optimize_method %in% c("ROI", "quadprog", "glpk", "symphony", "ipop", "cplex")){
obj_names <- unlist(lapply(portfolio$objectives, function(x) x$name))
if(any(obj_names %in% c("CVaR", "ES", "ETL"))){
.formals <- modify.args(formals=.formals, arglist=list(ROI=TRUE), dots=TRUE)
}
}
if("R" %in% names(.formals)) .formals <- modify.args(formals=.formals, arglist=NULL, R=R, dots=FALSE)
if("portfolio" %in% names(.formals)) .formals <- modify.args(formals=.formals, arglist=NULL, portfolio=portfolio, dots=FALSE)
.formals$... <- NULL
# call momentFUN
mout <- try(do.call(momentFUN, .formals), silent=TRUE)
if(inherits(mout, "try-error")) {
message(paste("portfolio moment function failed with message", mout))
} else {
#.args_env <- as.environment(mout)
#.args_env <- new.env()
# Assign each element of mout to the .args_env environment
#for(name in names(mout)){
# .args_env[[name]] <- mout[[name]]
#}
dotargs <- mout
}
# Function to normalize weights to min_sum and max_sum
# This function could be replaced by rp_transform
normalize_weights <- function(weights){
# normalize results if necessary
if(!is.null(constraints$min_sum) | !is.null(constraints$max_sum)){
# the user has passed in either min_sum or max_sum constraints for the portfolio, or both.
# we'll normalize the weights passed in to whichever boundary condition has been violated
# NOTE: this means that the weights produced by a numeric optimization algorithm like DEoptim
# might violate your constraints, so you'd need to renormalize them after optimizing
# we'll create functions for that so the user is less likely to mess it up.
##' NOTE: need to normalize in the optimization wrapper too before we return, since we've normalized in here
##' In Kris' original function, this was manifested as a full investment constraint
if(!is.null(constraints$max_sum) & constraints$max_sum != Inf ) {
max_sum=constraints$max_sum
if(sum(weights)>max_sum) { weights<-(max_sum/sum(weights))*weights } # normalize to max_sum
}
if(!is.null(constraints$min_sum) & constraints$min_sum != -Inf ) {
min_sum=constraints$min_sum
if(sum(weights)<min_sum) { weights<-(min_sum/sum(weights))*weights } # normalize to min_sum
}
} # end min_sum and max_sum normalization
return(weights)
}
# DEoptim optimization method
if(optimize_method == "DEoptim"){
stopifnot("package:DEoptim" %in% search() || require("DEoptim",quietly = TRUE))
# DEoptim does 200 generations by default, so lets set the size of each generation to search_size/200)
if(hasArg(itermax)) itermax=match.call(expand.dots=TRUE)$itermax else itermax=N*50
NP <- round(search_size/itermax)
if(NP < (N * 10)) NP <- N * 10
if(NP > 2000) NP <- 2000
if(!hasArg(itermax)) {
itermax <- round(search_size / NP)
if(itermax < 50) itermax <- 50 #set minimum number of generations
}
#check to see whether we need to disable foreach for parallel optimization, esp if called from inside foreach
if(hasArg(parallel)) parallel <- match.call(expand.dots=TRUE)$parallel else parallel <- TRUE
if(!isTRUE(parallel) && 'package:foreach' %in% search()){
registerDoSEQ()
}
DEcformals <- formals(DEoptim.control)
DEcargs <- names(DEcformals)
if( is.list(dotargs) ){
pm <- pmatch(names(dotargs), DEcargs, nomatch = 0L)
names(dotargs[pm > 0L]) <- DEcargs[pm]
DEcformals$NP <- NP
DEcformals$itermax <- itermax
DEcformals[pm] <- dotargs[pm > 0L]
if(!hasArg(strategy)) DEcformals$strategy=6 # use DE/current-to-p-best/1
if(!hasArg(reltol)) DEcformals$reltol=.000001 # 1/1000 of 1% change in objective is significant
if(!hasArg(steptol)) DEcformals$steptol=round(N*1.5) # number of assets times 1.5 tries to improve
if(!hasArg(c)) DEcformals$c=.4 # JADE mutation parameter, this could maybe use some adjustment
if(!hasArg(storepopfrom)) DEcformals$storepopfrom=1
if(isTRUE(parallel) && 'package:foreach' %in% search()){
if(!hasArg(parallelType) ) DEcformals$parallelType='auto' #use all cores
if(!hasArg(packages) ) DEcformals$packages <- names(sessionInfo()$otherPkgs) #use all packages
}
#TODO FIXME also check for a passed in controlDE list, including checking its class, and match formals
}
if(hasArg(traceDE)) traceDE=match.call(expand.dots=TRUE)$traceDE else traceDE=TRUE
DEcformals$trace <- traceDE
if(isTRUE(trace)) {
#we can't pass trace=TRUE into constrained objective with DEoptim, because it expects a single numeric return
tmptrace <- trace
assign('.objectivestorage', list(), pos='.GlobalEnv')
trace=FALSE
}
# get upper and lower weights parameters from constraints
upper <- constraints$max
lower <- constraints$min
# issue message if min_sum and max_sum are restrictive
if((constraints$max_sum - constraints$min_sum) < 0.02){
message("Leverage constraint min_sum and max_sum are restrictive,
consider relaxing. e.g. 'full_investment' constraint should be min_sum=0.99 and max_sum=1.01")
}
#if(hasArg(rpseed)){
# seed <- match.call(expand.dots=TRUE)$rpseed
# DEcformals$initialpop <- seed
# rpseed <- FALSE
#} else {
# rpseed <- TRUE
#}
#if(hasArg(rpseed) & isTRUE(rpseed)) {
# # initial seed population is generated with random_portfolios function
# # if(hasArg(eps)) eps=match.call(expand.dots=TRUE)$eps else eps = 0.01
# if(hasArg(rp_method)) rp_method=match.call(expand.dots=TRUE)$rp_method else rp_method="sample"
# if(hasArg(eliminate)) eliminate=match.call(expand.dots=TRUE)$eliminate else eliminate=TRUE
# if(hasArg(fev)) fev=match.call(expand.dots=TRUE)$fev else fev=0:5
# rp <- random_portfolios(portfolio=portfolio, permutations=NP, rp_method=rp_method, eliminate=eliminate, fev=fev)
# DEcformals$initialpop <- rp
#}
# Use rp as the initial population or generate from random portfolios
if(!is.null(rp)){
rp_len <- min(nrow(rp), NP)
seed <- rp[1:rp_len,]
DEcformals$initialpop <- seed
} else{
# Initial seed population is generated with random_portfolios function if rp is not passed in
if(hasArg(rp_method)) rp_method=match.call(expand.dots=TRUE)$rp_method else rp_method="sample"
# if(hasArg(eliminate)) eliminate=match.call(expand.dots=TRUE)$eliminate else eliminate=TRUE
if(hasArg(fev)) fev=match.call(expand.dots=TRUE)$fev else fev=0:5
rp <- random_portfolios(portfolio=portfolio, permutations=(NP+1), rp_method=rp_method, eliminate=FALSE, fev=fev)
DEcformals$initialpop <- rp
}
controlDE <- do.call(DEoptim.control, DEcformals)
# We are passing fn_map to the optional fnMap function to do the
# transformation so we need to force normalize=FALSE in call to constrained_objective
minw = try(DEoptim( constrained_objective, lower=lower[1:N], upper=upper[1:N], control=controlDE, R=R, portfolio=portfolio, env=dotargs, normalize=FALSE, fnMap=function(x) fn_map(x, portfolio=portfolio)$weights), silent=TRUE)
if(inherits(minw, "try-error")) {
message(minw)
minw=NULL
}
if(is.null(minw)){
message(paste("Optimizer was unable to find a solution for target"))
return (paste("Optimizer was unable to find a solution for target"))
}
if(isTRUE(tmptrace)) trace <- tmptrace
weights <- as.vector(minw$optim$bestmem)
print(weights)
# is it necessary to normalize the weights here?
# weights <- normalize_weights(weights)
names(weights) <- colnames(R)
obj_vals <- constrained_objective(w=weights, R=R, portfolio, trace=TRUE, normalize=FALSE, env=dotargs)$objective_measures
out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=minw$optim$bestval, call=call)
if (isTRUE(trace)){
out$DEoutput <- minw
out$DEoptim_objective_results <- try(get('.objectivestorage',pos='.GlobalEnv'),silent=TRUE)
rm('.objectivestorage',pos='.GlobalEnv')
}
} ## end case for DEoptim
# case for random portfolios optimization method
if(optimize_method=="random"){
# issue message if min_sum and max_sum are too restrictive
if((constraints$max_sum - constraints$min_sum) < 0.02){
message("Leverage constraint min_sum and max_sum are restrictive,
consider relaxing. e.g. 'full_investment' constraint should be min_sum=0.99 and max_sum=1.01")
}
#' call random_portfolios() with portfolio and search_size to create matrix of portfolios
if(missing(rp) | is.null(rp)){
if(hasArg(rp_method)) rp_method=match.call(expand.dots=TRUE)$rp_method else rp_method="sample"
if(hasArg(eliminate)) eliminate=match.call(expand.dots=TRUE)$eliminate else eliminate=TRUE
if(hasArg(fev)) fev=match.call(expand.dots=TRUE)$fev else fev=0:5
rp <- random_portfolios(portfolio=portfolio, permutations=search_size, rp_method=rp_method, eliminate=eliminate, fev=fev)
}
#' store matrix in out if trace=TRUE
if (isTRUE(trace)) out$random_portfolios <- rp
# rp is already being generated with a call to fn_map so set normalize=FALSE in the call to constrained_objective
#' write foreach loop to call constrained_objective() with each portfolio
if ("package:foreach" %in% search() & !hasArg(parallel)){
rp_objective_results <- foreach(ii=1:nrow(rp), .errorhandling='pass') %dopar% constrained_objective(w=rp[ii,], R=R, portfolio=portfolio, trace=trace, env=dotargs, normalize=FALSE)
} else {
rp_objective_results <- apply(rp, 1, constrained_objective, R=R, portfolio=portfolio, trace=trace, normalize=FALSE, env=dotargs)
}
#' if trace=TRUE , store results of foreach in out$random_results
if(isTRUE(trace)) out$random_portfolio_objective_results <- rp_objective_results
#' loop through results keeping track of the minimum value of rp_objective_results[[objective]]$out
search <- vector(length=length(rp_objective_results))
# first we construct the vector of results
for (i in 1:length(search)) {
if (isTRUE(trace)) {
search[i] <- ifelse(try(rp_objective_results[[i]]$out), rp_objective_results[[i]]$out,1e6)
} else {
search[i] <- as.numeric(rp_objective_results[[i]])
}
}
# now find the weights that correspond to the minimum score from the constrained objective
# and normalize_weights so that we meet our min_sum/max_sum constraints
# Is it necessary to normalize the weights at all with random portfolios?
if (isTRUE(trace)) {
min_objective_weights <- try(normalize_weights(rp_objective_results[[which.min(search)]]$weights))
} else {
min_objective_weights <- try(normalize_weights(rp[which.min(search),]))
}
#' re-call constrained_objective on the best portfolio, as above in DEoptim, with trace=TRUE to get results for out list
out$weights <- min_objective_weights
obj_vals <- try(constrained_objective(w=min_objective_weights, R=R, portfolio=portfolio, trace=TRUE, normalize=FALSE, env=dotargs)$objective_measures)
out$objective_measures <- obj_vals
out$opt_values <- obj_vals
out$call <- call
#' construct out list to be as similar as possible to DEoptim list, within reason
} ## end case for random
roi_solvers <- c("ROI", "quadprog", "glpk", "symphony", "ipop", "cplex")
if(optimize_method %in% roi_solvers){
# check for a control argument for list of solver control arguments
if(hasArg(control)) control=match.call(expand.dots=TRUE)$control else control=NULL
# This takes in a regular portfolio object and extracts the constraints and
# objectives and converts them for input. to a closed form solver using
# ROI as an interface.
moments <- list(mean=rep(0, N))
alpha <- 0.05
if(!is.null(constraints$return_target)){
target <- constraints$return_target
} else {
target <- NA
}
lambda <- 1
# list of valid objective names for ROI solvers
valid_objnames <- c("HHI", "mean", "var", "sd", "StdDev", "CVaR", "ES", "ETL")
#objnames <- unlist(lapply(portfolio$objectives, function(x) x$name))
for(objective in portfolio$objectives){
if(objective$enabled){
if(!(objective$name %in% valid_objnames)){
stop("ROI only solves mean, var/StdDev, HHI, or sample ETL/ES/CVaR type business objectives, choose a different optimize_method.")
}
# Grab the arguments list per objective
# Currently we are only getting arguments for "p" and "clean", not sure if we need others for the ROI QP/LP solvers
# if(length(objective$arguments) >= 1) arguments <- objective$arguments else arguments <- list()
arguments <- objective$arguments
if(!is.null(arguments$clean)) clean <- arguments$clean else clean <- "none"
# Note: arguments$p grabs arguments$portfolio_method if no p is specified
# so we need to be explicit with arguments[["p"]]
if(!is.null(arguments[["p"]])) alpha <- arguments$p else alpha <- alpha
if(alpha > 0.5) alpha <- (1 - alpha)
# Some of the sub-functions for optimizations use the returns object as
# part of the constraints matrix (e.g. etl_opt and etl_milp_opt) so we
# will store the cleaned returns in the moments object. This may not
# be the most efficient way to pass around a cleaned returns object,
# but it will keep it separate from the R object passed in by the user
# and avoid "re-cleaning" already cleaned returns if specified in
# multiple objectives.
if(clean != "none") moments$cleanR <- Return.clean(R=R, method=clean)
# Use $mu and $sigma estimates from momentFUN if available, fall back to
# calculating sample mean and variance
if(objective$name == "mean"){
if(!is.null(mout$mu)){
moments[["mean"]] <- as.vector(mout$mu)
} else {
moments[["mean"]] <- try(as.vector(apply(Return.clean(R=R, method=clean), 2, "mean", na.rm=TRUE)), silent=TRUE)
}
} else if(objective$name %in% c("StdDev", "sd", "var")){
if(!is.null(mout$sigma)){
moments[["var"]] <- mout$sigma
} else {
moments[["var"]] <- try(var(x=Return.clean(R=R, method=clean), na.rm=TRUE), silent=TRUE)
}
} else if(objective$name %in% c("CVaR", "ES", "ETL")){
# do nothing (no point in computing ES here)
moments[[objective$name]] <- ""
} else {
# I'm not sure this is serving any purpose since we control the types
# of problems solved with ROI. The min ES problem only uses
# moments$mu if a target return is specified.
moments[[objective$name]] <- try(eval(as.symbol(objective$name))(Return.clean(R=R, method=clean)), silent=TRUE)
}
target <- ifelse(!is.null(objective$target), objective$target, target)
# alpha <- ifelse(!is.null(objective$alpha), objective$alpha, alpha)
# only accept confidence level for ES/ETL/CVaR to come from the
# arguments list to be consistent with how this is done in other solvers.
lambda <- ifelse(!is.null(objective$risk_aversion), objective$risk_aversion, lambda)
if(!is.null(objective$conc_aversion)) lambda_hhi <- objective$conc_aversion else lambda_hhi <- NULL
if(!is.null(objective$conc_groups)) conc_groups <- objective$conc_groups else conc_groups <- NULL
}
}
if("var" %in% names(moments)){
# Set a default solver if optimize_method == "ROI", otherwise assume the
# optimize_method specified by the user is the solver for ROI
if(optimize_method == "ROI"){
solver <- "quadprog"
} else {
solver <- optimize_method
}
# Minimize variance if the only objective specified is variance
# Maximize Quadratic Utility if var and mean are specified as objectives
if(!is.null(constraints$turnover_target) | !is.null(constraints$ptc)){
if(!is.null(constraints$turnover_target) & !is.null(constraints$ptc)){
warning("Turnover and proportional transaction cost constraints detected, only running optimization for turnover constraint.")
constraints$ptc <- NULL
}
if(!is.null(constraints$turnover_target) & is.null(constraints$ptc)){
qp_result <- gmv_opt_toc(R=R, constraints=constraints, moments=moments, lambda=lambda, target=target, init_weights=portfolio$assets, solver=solver, control=control)
weights <- qp_result$weights
# obj_vals <- constrained_objective(w=weights, R=R, portfolio, trace=TRUE, normalize=FALSE)$objective_measures
obj_vals <- qp_result$obj_vals
out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=qp_result$out, call=call)
}
if(!is.null(constraints$ptc) & is.null(constraints$turnover_target)){
qp_result <- gmv_opt_ptc(R=R, constraints=constraints, moments=moments, lambda=lambda, target=target, init_weights=portfolio$assets, solver=solver, control=control)
weights <- qp_result$weights
# obj_vals <- constrained_objective(w=weights, R=R, portfolio, trace=TRUE, normalize=FALSE)$objective_measures
obj_vals <- qp_result$obj_vals
out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=qp_result$out, call=call)
}
} else {
# if(hasArg(ef)) ef=match.call(expand.dots=TRUE)$ef else ef=FALSE
if(hasArg(maxSR)) maxSR=match.call(expand.dots=TRUE)$maxSR else maxSR=FALSE
if(maxSR){
target <- max_sr_opt(R=R, constraints=constraints, moments=moments, lambda_hhi=lambda_hhi, conc_groups=conc_groups, solver=solver, control=control)
# need to set moments$mean=0 here because quadratic utility and target return is sensitive to returning no solution
tmp_moments_mean <- moments$mean
moments$mean <- rep(0, length(moments$mean))
}
roi_result <- gmv_opt(R=R, constraints=constraints, moments=moments, lambda=lambda, target=target, lambda_hhi=lambda_hhi, conc_groups=conc_groups, solver=solver, control=control)
weights <- roi_result$weights
# obj_vals <- constrained_objective(w=weights, R=R, portfolio, trace=TRUE, normalize=FALSE)$objective_measures
obj_vals <- roi_result$obj_vals
if(maxSR){
# need to recalculate mean here if we are maximizing sharpe ratio
port.mean <- as.numeric(sum(weights * tmp_moments_mean))
names(port.mean) <- "mean"
obj_vals$mean <- port.mean
}
out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=roi_result$out, call=call)
}
}
if(length(names(moments)) == 1 & "mean" %in% names(moments)) {
# Set a default solver if optimize_method == "ROI", otherwise assume the
# optimize_method specified by the user is the solver for ROI
if(optimize_method == "ROI"){
solver <- "glpk"
} else {
solver <- optimize_method
}
# Maximize return if the only objective specified is mean
if(!is.null(constraints$max_pos) | !is.null(constraints$leverage)) {
# This is an MILP problem if max_pos is specified as a constraint
roi_result <- maxret_milp_opt(R=R, constraints=constraints, moments=moments, target=target, solver=solver, control=control)
weights <- roi_result$weights
# obj_vals <- constrained_objective(w=weights, R=R, portfolio, trace=TRUE, normalize=FALSE)$objective_measures
obj_vals <- roi_result$obj_vals
out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=roi_result$out, call=call)
} else {
# Maximize return LP problem
roi_result <- maxret_opt(R=R, constraints=constraints, moments=moments, target=target, solver=solver, control=control)
weights <- roi_result$weights
# obj_vals <- constrained_objective(w=weights, R=R, portfolio, trace=TRUE, normalize=FALSE)$objective_measures
obj_vals <- roi_result$obj_vals
out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=roi_result$out, call=call)
}
}
if( any(c("CVaR", "ES", "ETL") %in% names(moments)) ) {
# Set a default solver if optimize_method == "ROI", otherwise assume the
# optimize_method specified by the user is the solver for ROI
if(optimize_method == "ROI"){
solver <- "glpk"
} else {
solver <- optimize_method
}
if(hasArg(ef)) ef=match.call(expand.dots=TRUE)$ef else ef=FALSE
if(hasArg(maxSTARR)) maxSTARR=match.call(expand.dots=TRUE)$maxSTARR else maxSTARR=TRUE
if(ef) meanetl <- TRUE else meanetl <- FALSE
tmpnames <- c("CVaR", "ES", "ETL")
idx <- which(tmpnames %in% names(moments))
# Minimize sample ETL/ES/CVaR if CVaR, ETL, or ES is specified as an objective
if(length(moments) == 2 & all(moments$mean != 0) & ef==FALSE & maxSTARR){
# This is called by meanetl.efficient.frontier and we do not want that for efficient frontiers, need to have ef==FALSE
target <- mean_etl_opt(R=R, constraints=constraints, moments=moments, alpha=alpha, solver=solver, control=control)
meanetl <- TRUE
}
if(!is.null(constraints$max_pos)) {
# This is an MILP problem if max_pos is specified as a constraint
roi_result <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=target, alpha=alpha, solver=solver, control=control)
weights <- roi_result$weights
# obj_vals <- constrained_objective(w=weights, R=R, portfolio, trace=TRUE, normalize=FALSE)$objective_measures
# obj_vals <- roi_result$obj_vals
# calculate obj_vals based on solver output
obj_vals <- list()
if(meanetl) obj_vals$mean <- sum(weights * moments$mean)
obj_vals[[tmpnames[idx]]] <- roi_result$out
out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=roi_result$out, call=call)
} else {
# Minimize sample ETL/ES/CVaR LP Problem
roi_result <- etl_opt(R=R, constraints=constraints, moments=moments, target=target, alpha=alpha, solver=solver, control=control)
weights <- roi_result$weights
# obj_vals <- constrained_objective(w=weights, R=R, portfolio, trace=TRUE, normalize=FALSE)$objective_measures
# obj_vals <- roi_result$obj_vals
obj_vals <- list()
if(meanetl) obj_vals$mean <- sum(weights * moments$mean)
obj_vals[[tmpnames[idx]]] <- roi_result$out
out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=roi_result$out, call=call)
}
}
# Set here at the end so we get optimize.portfolio.ROI and not optimize.portfolio.{solver} classes
optimize_method <- "ROI"
} ## end case for ROI
## case if method=pso---particle swarm
if(optimize_method=="pso"){
stopifnot("package:pso" %in% search() || require("pso",quietly = TRUE) )
if(hasArg(maxit)) maxit=match.call(expand.dots=TRUE)$maxit else maxit=N*50
controlPSO <- list(trace=FALSE, fnscale=1, maxit=1000, maxf=Inf, abstol=-Inf, reltol=0)
PSOcargs <- names(controlPSO)
if( is.list(dotargs) ){
pm <- pmatch(names(dotargs), PSOcargs, nomatch = 0L)
names(dotargs[pm > 0L]) <- PSOcargs[pm]
controlPSO$maxit <- maxit
controlPSO[pm] <- dotargs[pm > 0L]
if(!hasArg(reltol)) controlPSO$reltol <- .000001 # 1/1000 of 1% change in objective is significant
if(hasArg(trace) && try(trace==TRUE,silent=TRUE)) controlPSO$trace <- TRUE
if(hasArg(trace) && isTRUE(trace)) {
controlPSO$trace <- TRUE
controlPSO$trace.stats=TRUE
}
}
# get upper and lower weights parameters from constraints
upper <- constraints$max
lower <- constraints$min
minw <- try(psoptim( par = rep(NA, N), fn = constrained_objective, R=R, portfolio=portfolio, env=dotargs,
lower = lower[1:N] , upper = upper[1:N] , control = controlPSO)) # add ,silent=TRUE here?
if(inherits(minw,"try-error")) { minw=NULL }
if(is.null(minw)){
message(paste("Optimizer was unable to find a solution for target"))
return (paste("Optimizer was unable to find a solution for target"))
}
weights <- as.vector( minw$par)
weights <- normalize_weights(weights)
names(weights) <- colnames(R)
obj_vals <- constrained_objective(w=weights, R=R, portfolio=portfolio, trace=TRUE, env=dotargs)$objective_measures
out <- list(weights=weights,
objective_measures=obj_vals,
opt_values=obj_vals,
out=minw$value,
call=call)
if (isTRUE(trace)){
out$PSOoutput=minw
}
} ## end case for pso
## case if method=GenSA---Generalized Simulated Annealing
if(optimize_method=="GenSA"){
stopifnot("package:GenSA" %in% search() || require("GenSA",quietly = TRUE) )
if(hasArg(maxit)) maxit=match.call(expand.dots=TRUE)$maxit else maxit=N*50
controlGenSA <- list(maxit = 5000, threshold.stop = NULL, temp = 5230,
visiting.param = 2.62, acceptance.param = -5, max.time = NULL,
nb.stop.improvement = 1e+06, smooth = TRUE, max.call = 1e+07,
verbose = FALSE)
GenSAcargs <- names(controlGenSA)
if( is.list(dotargs) ){
pm <- pmatch(names(dotargs), GenSAcargs, nomatch = 0L)
names(dotargs[pm > 0L]) <- GenSAcargs[pm]
controlGenSA$maxit <- maxit
controlGenSA[pm] <- dotargs[pm > 0L]
if(hasArg(trace) && try(trace==TRUE,silent=TRUE)) controlGenSA$verbose <- TRUE
}
upper <- constraints$max
lower <- constraints$min
minw = try(GenSA( par = rep(1/N, N), lower = lower[1:N] , upper = upper[1:N], control = controlGenSA,
fn = constrained_objective , R=R, portfolio=portfolio, env=dotargs)) # add ,silent=TRUE here?
if(inherits(minw,"try-error")) { minw=NULL }
if(is.null(minw)){
message(paste("Optimizer was unable to find a solution for target"))
return (paste("Optimizer was unable to find a solution for target"))
}
weights <- as.vector(minw$par)
weights <- normalize_weights(weights)
names(weights) <- colnames(R)
obj_vals <- constrained_objective(w=weights, R=R, portfolio=portfolio, trace=TRUE, env=dotargs)$objective_measures
out = list(weights=weights,
objective_measures=obj_vals,
opt_values=obj_vals,
out=minw$value,
call=call)
if (isTRUE(trace)){
out$GenSAoutput=minw
}
} ## end case for GenSA
# Prepare for final object to return
end_t <- Sys.time()
# print(c("elapsed time:",round(end_t-start_t,2),":diff:",round(diff,2), ":stats: ", round(out$stats,4), ":targets:",out$targets))
if(message) message(c("elapsed time:", end_t-start_t))
out$portfolio <- portfolio
if(trace) out$R <- R
out$data_summary <- list(first=first(R), last=last(R))
out$elapsed_time <- end_t - start_t
out$end_t <- as.character(Sys.time())
# return a $regime element to indicate what regime portfolio used for
# optimize.portfolio. The regime information is used in extractStats and
# extractObjectiveMeasures
if(regime.switching){
out$regime <- regime.idx
}
class(out) <- c(paste("optimize.portfolio", optimize_method, sep='.'), "optimize.portfolio")
return(out)
}
#' Constrained optimization of portfolios
#'
#' This function aims to provide a wrapper for constrained optimization of
#' portfolios that specify constraints and objectives.
#'
#' @details
#' This function currently supports DEoptim, random portfolios, pso, GenSA, and ROI as back ends.
#' Additional back end contributions for Rmetrics, ghyp, etc. would be welcome.
#'
#' When using random portfolios, search_size is precisely that, how many
#' portfolios to test. You need to make sure to set your feasible weights
#' in generatesequence to make sure you have search_size unique
#' portfolios to test, typically by manipulating the 'by' parameter
#' to select something smaller than .01
#' (I often use .002, as .001 seems like overkill)
#'
#' When using DE, search_size is decomposed into two other parameters
#' which it interacts with, NP and itermax.
#'
#' NP, the number of members in each population, is set to cap at 2000 in
#' DEoptim, and by default is the number of parameters (assets/weights) * 10.
#'
#' itermax, if not passed in dots, defaults to the number of parameters (assets/weights) * 50.
#'
#' When using GenSA and want to set \code{verbose=TRUE}, instead use \code{trace}.
#'
#' If \code{optimize_method="ROI"} is specified, a default solver will be
#' selected based on the optimization problem. The \code{glpk} solver is the
#' default solver for LP and MILP optimization problems. The \code{quadprog}
#' solver is the default solver for QP optimization problems. For example,
#' \code{optimize_method = "quadprog"} can be specified and the optimization
#' problem will be solved via ROI using the quadprog solver.
#'
#' The extension to ROI solves a limited type of convex optimization problems:
#' \itemize{
#' \item{Maxmimize portfolio return subject leverage, box, group, position limit, target mean return, and/or factor exposure constraints on weights.}
#' \item{Minimize portfolio variance subject to leverage, box, group, turnover, and/or factor exposure constraints (otherwise known as global minimum variance portfolio).}
#' \item{Minimize portfolio variance subject to leverage, box, group, and/or factor exposure constraints and a desired portfolio return.}
#' \item{Maximize quadratic utility subject to leverage, box, group, target mean return, turnover, and/or factor exposure constraints and risk aversion parameter.
#' (The risk aversion parameter is passed into \code{optimize.portfolio} as an added argument to the \code{portfolio} object).}
#' \item{Maximize portfolio mean return per unit standard deviation (i.e. the Sharpe Ratio) can be done by specifying \code{maxSR=TRUE} in \code{optimize.portfolio}.
#' If both mean and StdDev are specified as objective names, the default action is to maximize quadratic utility, therefore \code{maxSR=TRUE} must be specified to maximize Sharpe Ratio.}
#' \item{Minimize portfolio ES/ETL/CVaR optimization subject to leverage, box, group, position limit, target mean return, and/or factor exposure constraints and target portfolio return.}
#' \item{Maximize portfolio mean return per unit ES/ETL/CVaR (i.e. the STARR Ratio) can be done by specifying \code{maxSTARR=TRUE} in \code{optimize.portfolio}.
#' If both mean and ES/ETL/CVaR are specified as objective names, the default action is to maximize mean return per unit ES/ETL/CVaR.}
#' }
#' These problems also support a weight_concentration objective where concentration
#' of weights as measured by HHI is added as a penalty term to the quadratic objective.
#'
#' Because these convex optimization problem are standardized, there is no need for a penalty term.
#' The \code{multiplier} argument in \code{\link{add.objective}} passed into the complete constraint object are ingnored by the ROI solver.
#'
#' @note
#' An object of class \code{v1_constraint} can be passed in for the \code{constraints} argument.
#' The \code{v1_constraint} object was used in the previous 'v1' specification to specify the
#' constraints and objectives for the optimization problem, see \code{\link{constraint}}.
#' We will attempt to detect if the object passed into the constraints argument
#' is a \code{v1_constraint} object and update to the 'v2' specification by adding the
#' constraints and objectives to the \code{portfolio} object.
#'
#' @param R an xts, vector, matrix, data frame, timeSeries or zoo object of asset returns
#' @param portfolio an object of type "portfolio" specifying the constraints and objectives for the optimization
#' @param constraints default=NULL, a list of constraint objects. An object of class 'v1_constraint' can be passed in here.
#' @param objectives default=NULL, a list of objective objects.
#' @param optimize_method one of "DEoptim", "random", "ROI", "pso", "GenSA". A solver
#' for ROI can also be specified and will be solved using ROI. See Details.
#' @param search_size integer, how many portfolios to test, default 20,000
#' @param trace TRUE/FALSE if TRUE will attempt to return additional information on the path or portfolios searched
#' @param \dots any other passthru parameters
#' @param rp matrix of random portfolio weights, default NULL, mostly for automated use by rebalancing optimization or repeated tests on same portfolios
#' @param momentFUN the name of a function to call to set portfolio moments, default \code{\link{set.portfolio.moments_v2}}
#' @param message TRUE/FALSE. The default is message=FALSE. Display messages if TRUE.
#'
#' @return a list containing the following elements
#' \itemize{
#' \item{\code{weights}:}{ The optimal set weights.}
#' \item{\code{objective_measures}:}{ A list containing the value of each objective corresponding to the optimal weights.}
#' \item{\code{opt_values}:}{ A list containing the value of each objective corresponding to the optimal weights.}
#' \item{\code{out}:}{ The output of the solver.}
#' \item{\code{call}:}{ The function call.}
#' \item{\code{portfolio}:}{ The portfolio object.}
#' \item{\code{R}:}{ The asset returns.}
#' \item{\code{data summary:}}{ The first row and last row of \code{R}.}
#' \item{\code{elapsed_time:}}{ The amount of time that elapses while the optimization is run.}
#' \item{\code{end_t:}}{ The date and time the optimization completed.}
#' }
#' When Trace=TRUE is specified, the following elements will be returned in
#' addition to the elements above. The output depends on the optimization
#' method and is specific to each solver. Refer to the documentation of the
#' desired solver for more information.
#'
#' \code{optimize_method="random"}
#' \itemize{
#' \item{\code{random_portfolios}:}{ A matrix of the random portfolios.}
#' \item{\code{random_portfolio_objective_results}:}{ A list of the following elements for each random portfolio.}
#' \itemize{
#' \item{\code{out}:}{ The output value of the solver corresponding to the random portfolio weights.}
#' \item{\code{weights}:}{ The weights of the random portfolio.}
#' \item{\code{objective_measures}:}{ A list of each objective measure corresponding to the random portfolio weights.}
#' }
#' }
#'
#' \code{optimize_method="DEoptim"}
#' \itemize{
#' \item{\code{DEoutput:}}{ A list (of length 2) containing the following elements:}
#' \itemize{
#' \item{\code{optim}}
#' \item{\code{member}}
#' }
#' \item{\code{DEoptim_objective_results}:}{ A list containing the following elements for each intermediate population.}
#' \itemize{
#' \item{\code{out}:}{ The output of the solver.}
#' \item{\code{weights}:}{ Population weights.}
#' \item{\code{init_weights}:}{ Initial population weights.}
#' \item{\code{objective_measures}:}{ A list of each objective measure corresponding to the weights}
#' }
#' }
#'
#' \code{optimize_method="pso"}
#' \itemize{
#' \item{\code{PSOoutput}:}{ A list containing the following elements:}
#' \itemize{
#' \item{par}
#' \item{value}
#' \item{counts}
#' \item{convergence}
#' \item{message}
#' \item{stats}
#' }
#' }
#'
#' \code{optimize_method="GenSA"}
#' \itemize{
#' \item{\code{GenSAoutput:}}{ A list containing the following elements:}
#' \itemize{
#' \item{value}
#' \item{par}
#' \item{trace.mat}
#' \item{counts}
#' }
#' }
#'
#' @author Kris Boudt, Peter Carl, Brian G. Peterson, Ross Bennett
#' @aliases optimize.portfolio_v2 optimize.portfolio_v1
#' @seealso \code{\link{portfolio.spec}}
#' @name optimize.portfolio
#' @export
optimize.portfolio <- optimize.portfolio_v2
#' @rdname optimize.portfolio.rebalancing
#' @name optimize.portfolio.rebalancing
#' @export
optimize.portfolio.rebalancing_v1 <- function(R,constraints,optimize_method=c("DEoptim","random","ROI"), search_size=20000, trace=FALSE, ..., rp=NULL, rebalance_on=NULL, training_period=NULL, trailing_periods=NULL)
{
stopifnot("package:foreach" %in% search() || require("foreach",quietly=TRUE))
start_t<-Sys.time()
#store the call for later
call <- match.call()
if(optimize_method=="random"){
#' call random_portfolios() with constraints and search_size to create matrix of portfolios
if(is.null(rp))
rp<-random_portfolios(rpconstraints=constraints,permutations=search_size)
} else {
rp=NULL
}
if(is.null(training_period)) {if(nrow(R)<36) training_period=nrow(R) else training_period=36}
if (is.null(trailing_periods)){
# define the index endpoints of our periods
ep.i<-endpoints(R,on=rebalance_on)[which(endpoints(R, on = rebalance_on)>=training_period)]
# now apply optimize.portfolio to the periods, in parallel if available
out_list<-foreach(ep=iter(ep.i), .errorhandling='pass', .packages='PortfolioAnalytics') %dopar% {
optimize.portfolio(R[1:ep,],constraints=constraints,optimize_method=optimize_method, search_size=search_size, trace=trace, rp=rp, parallel=FALSE, ...=...)
}
} else {
# define the index endpoints of our periods
ep.i<-endpoints(R,on=rebalance_on)[which(endpoints(R, on = rebalance_on)>=training_period)]
# now apply optimize.portfolio to the periods, in parallel if available
out_list<-foreach(ep=iter(ep.i), .errorhandling='pass', .packages='PortfolioAnalytics') %dopar% {
optimize.portfolio(R[(ifelse(ep-trailing_periods>=1,ep-trailing_periods,1)):ep,],constraints=constraints,optimize_method=optimize_method, search_size=search_size, trace=trace, rp=rp, parallel=FALSE, ...=...)
}
}
names(out_list)<-index(R[ep.i])
end_t<-Sys.time()
message(c("overall elapsed time:",end_t-start_t))
class(out_list)<-c("optimize.portfolio.rebalancing")
return(out_list)
}
#' Portfolio Optimization with Rebalancing Periods
#'
#' Portfolio optimization with support for rebalancing periods for
#' out-of-sample testing (i.e. backtesting)
#'
#' @details
#' Run portfolio optimization with periodic rebalancing at specified time periods.
#' Running the portfolio optimization with periodic rebalancing can help
#' refine the constraints and objectives by evaluating the out of sample
#' performance of the portfolio based on historical data
#'
#' This function is a essentially a wrapper around \code{optimize.portfolio}
#' and thus the discussion in the Details section of the
#' \code{optimize.portfolio} help file is valid here as well.
#'
#' This function is massively parallel and requires the 'foreach' package. It
#' is suggested to register a parallel backend.
#'
#' @param R an xts, vector, matrix, data frame, timeSeries or zoo object of asset returns
#' @param portfolio an object of type "portfolio" specifying the constraints
#' and objectives for the optimization
#' @param constraints default NULL, a list of constraint objects
#' @param objectives default NULL, a list of objective objects
#' @param optimize_method one of "DEoptim", "random", "pso", "GenSA", or "ROI"
#' @param search_size integer, how many portfolios to test, default 20,000
#' @param trace TRUE/FALSE if TRUE will attempt to return additional
#' information on the path or portfolios searched
#' @param \dots any other passthru parameters to \code{\link{optimize.portfolio}}
#' @param rp a set of random portfolios passed into the function to prevent recalculation
#' @param rebalance_on character string of period to rebalance on. See
#' \code{\link[xts]{endpoints}} for valid names.
#' @param training_period an integer of the number of periods to use as
#' a training data in the front of the returns data
#' @param trailing_periods an integer with the number of periods to roll over
#' (i.e. width of the moving or rolling window), the default is NULL will
#' run using the returns data from inception
#' @return a list containing the following elements
#' \itemize{
#' \item{\code{portfolio}:}{ The portfolio object.}
#' \item{\code{R}:}{ The asset returns.}
#' \item{\code{call}:}{ The function call.}
#' \item{\code{elapsed_time:}}{ The amount of time that elapses while the
#' optimization is run.}
#' \item{\code{opt_rebalancing:}}{ A list of \code{optimize.portfolio}
#' objects computed at each rebalancing period.}
#' }
#' @author Kris Boudt, Peter Carl, Brian G. Peterson
#' @name optimize.portfolio.rebalancing
#' @aliases optimize.portfolio.rebalancing optimize.portfolio.rebalancing_v1
#' @seealso \code{\link{portfolio.spec}} \code{\link{optimize.portfolio}}
#' @examples
#' \dontrun{
#' data(edhec)
#' R <- edhec[,1:4]
#' funds <- colnames(R)
#'
#' portf <- portfolio.spec(funds)
#' portf <- add.constraint(portf, type="full_investment")
#' portf <- add.constraint(portf, type="long_only")
#' portf <- add.objective(portf, type="risk", name="StdDev")
#'
#' # Quarterly rebalancing with 5 year training period
#' bt.opt1 <- optimize.portfolio.rebalancing(R, portf,
#' optimize_method="ROI",
#' rebalance_on="quarters",
#' training_period=60)
#'
#' # Monthly rebalancing with 5 year training period and 4 year trailing (moving window)
#' bt.opt2 <- optimize.portfolio.rebalancing(R, portf,
#' optimize_method="ROI",
#' rebalance_on="months",
#' training_period=60,
#' trailing_period=48)
#' }
#' @export
optimize.portfolio.rebalancing <- function(R, portfolio=NULL, constraints=NULL, objectives=NULL, optimize_method=c("DEoptim","random","ROI"), search_size=20000, trace=FALSE, ..., rp=NULL, rebalance_on=NULL, training_period=NULL, trailing_periods=NULL)
{
stopifnot("package:foreach" %in% search() || require("foreach",quietly=TRUE))
stopifnot("package:iterators" %in% search() || require("iterators",quietly=TRUE))
# This is the case where the user has passed in a list of portfolio objects
# for the portfolio argument.
# Loop through the portfolio list and recursively call
# optimize.portfolio.rebalancing.
#Note that I return at the end of this block. I know it is not good practice
# to return before the end of a function, but I am not sure of another way
# to handle a list of portfolio objects with the recursive call to
# optimize.portfolio.
if(inherits(portfolio, "portfolio.list")){
n.portf <- length(portfolio)
opt.list <- vector("list", n.portf)
for(i in 1:length(opt.list)){
if(hasArg(message)) message=match.call(expand.dots=TRUE)$message else message=FALSE
if(message) cat("Starting optimization of portfolio ", i, "\n")
opt.list[[i]] <- optimize.portfolio.rebalancing(R=R,
portfolio=portfolio[[i]],
constraints=constraints,
objectives=objectives,
optimize_method=optimize_method,
search_size=search_size,
trace=trace,
...=...,
rp=rp,
rebalance_on=rebalance_on,
training_period=training_period,
trailing_periods=trailing_periods)
}
out <- combine.optimizations(opt.list)
class(out) <- "opt.rebal.list"
##### return here for portfolio.list because this is a recursive call
##### for optimize.portfolio.rebalancing
return(out)
}
# Store the call to return later
call <- match.call()
start_t<-Sys.time()
if (!is.null(portfolio) & !is.portfolio(portfolio)){
stop("you must pass in an object of class 'portfolio' to control the optimization")
}
if(hasArg(message)) message=match.call(expand.dots=TRUE)$message else message=FALSE
# Check for constraints and objectives passed in separately outside of the portfolio object
if(!is.null(constraints)){
if(inherits(constraints, "v1_constraint")){
if(is.null(portfolio)){
# If the user has not passed in a portfolio, we will create one for them
tmp_portf <- portfolio.spec(assets=constraints$assets)
}
message("constraint object passed in is a 'v1_constraint' object, updating to v2 specification")
portfolio <- update_constraint_v1tov2(portfolio=tmp_portf, v1_constraint=constraints)
# print.default(portfolio)
}
if(!inherits(constraints, "v1_constraint")){
# Insert the constraints into the portfolio object
portfolio <- insert_constraints(portfolio=portfolio, constraints=constraints)
}
}
if(!is.null(objectives)){
# Insert the objectives into the portfolio object
portfolio <- insert_objectives(portfolio=portfolio, objectives=objectives)
}
#store the call for later
call <- match.call()
if(optimize_method=="random"){
# get any rp related arguments passed in through dots
if(hasArg(rp_method)) rp_method=match.call(expand.dots=TRUE)$rp_method else rp_method="sample"
if(hasArg(eliminate)) eliminate=match.call(expand.dots=TRUE)$eliminate else eliminate=TRUE
if(hasArg(fev)) fev=match.call(expand.dots=TRUE)$fev else fev=0:5
#' call random_portfolios() with constraints and search_size to create matrix of portfolios
if(is.null(rp))
if(inherits(portfolio, "regime.portfolios")){
rp <- rp.regime.portfolios(regime=portfolio, permutations=search_size, rp_method=rp_method, eliminate=eliminate, fev=fev)
} else {
rp <- random_portfolios(portfolio=portfolio, permutations=search_size, rp_method=rp_method, eliminate=eliminate, fev=fev)
}
} else {
rp = NULL
}
if(is.null(training_period)) {if(nrow(R)<36) training_period=nrow(R) else training_period=36}
if (is.null(trailing_periods)){
# define the index endpoints of our periods
ep.i<-endpoints(R,on=rebalance_on)[which(endpoints(R, on = rebalance_on)>=training_period)]
# now apply optimize.portfolio to the periods, in parallel if available
out_list<-foreach(ep=iter(ep.i), .errorhandling='pass', .packages='PortfolioAnalytics') %dopar% {
optimize.portfolio(R[1:ep,], portfolio=portfolio, optimize_method=optimize_method, search_size=search_size, trace=trace, rp=rp, parallel=FALSE, ...=...)
}
} else {
# define the index endpoints of our periods
ep.i<-endpoints(R,on=rebalance_on)[which(endpoints(R, on = rebalance_on)>=training_period)]
# now apply optimize.portfolio to the periods, in parallel if available
out_list<-foreach(ep=iter(ep.i), .errorhandling='pass', .packages='PortfolioAnalytics') %dopar% {
optimize.portfolio(R[(ifelse(ep-trailing_periods>=1,ep-trailing_periods,1)):ep,], portfolio=portfolio, optimize_method=optimize_method, search_size=search_size, trace=trace, rp=rp, parallel=FALSE, ...=...)
}
}
# out_list is a list where each element is an optimize.portfolio object
# at each rebalance date
names(out_list)<-index(R[ep.i])
end_t <- Sys.time()
elapsed_time <- end_t - start_t
if(message) message(c("overall elapsed time:", end_t-start_t))
# out object to return
out <- list()
out$portfolio <- portfolio
out$R <- R
out$call <- call
out$elapsed_time <- elapsed_time
out$opt_rebalancing <- out_list
class(out) <- c("optimize.portfolio.rebalancing")
return(out)
}
#'execute multiple optimize.portfolio calls, presumably in parallel
#'
#' TODO write function to check sensitivity of optimal results by using optimize.portfolio.parallel results
#'
#' This function will not speed up optimization!
#'
#' This function exists to run multiple copies of optimize.portfolio, presumabley in parallel using foreach.
#'
#' This is typically done to test your parameter settings, specifically
#' total population size, but also possibly to help tune your
#' convergence settings, number of generations, stopping criteria,
#' etc.
#'
#' If you want to use all the cores on your multi-core computer, use
#' the parallel version of the apppropriate optimization engine, not
#' this function.
#'
#' @param R an xts, vector, matrix, data frame, timeSeries or zoo object of asset returns
#' @param constraints an object of type "constraints" specifying the constraints for the optimization, see \code{\link{constraint}}
#' @param optimize_method one of "DEoptim" or "random"
#' @param search_size integer, how many portfolios to test, default 20,000
#' @param trace TRUE/FALSE if TRUE will attempt to return additional information on the path or portfolios searched
#' @param \dots any other passthru parameters
#' @param nodes how many processes to run in the foreach loop, default 4
#'
#' @return a list containing the optimal weights, some summary statistics, the function call, and optionally trace information
#' @author Kris Boudt, Peter Carl, Brian G. Peterson
#' @export
optimize.portfolio.parallel <- function(R,constraints,optimize_method=c("DEoptim","random"), search_size=20000, trace=FALSE, ..., nodes=4)
{
stopifnot("package:foreach" %in% search() || require("foreach",quietly=TRUE))
optimize_method=optimize_method[1]
start_t<-Sys.time()
#store the call for later
call <- match.call()
opt_out_list<-foreach(1:nodes, packages='PortfolioAnalytics') %dopar% optimize.portfolio(R=R,constraints=constraints,optimize_method=optimize_method, search_size=search_size, trace=trace, ...)
end_t<-Sys.time()
message(c("overall elapsed time:",end_t-start_t))
class(opt_out_list)<-c("optimize.portfolio.parallel")
return(opt_out_list)
}
#TODO write function to compute an efficient frontier of optimal portfolios
###############################################################################
# $Id$
###############################################################################
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.