[Returnanalytics-commits] r2538 - pkg/PortfolioAnalytics/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jul 11 05:40:14 CEST 2013
Author: rossbennett34
Date: 2013-07-11 05:40:14 +0200 (Thu, 11 Jul 2013)
New Revision: 2538
Modified:
pkg/PortfolioAnalytics/R/constrained_objective.R
pkg/PortfolioAnalytics/R/optimize.portfolio.R
Log:
Adding constrained_objective_v2 to accept a portfolio object. Modified optimize.portfolio_v2 for DEoptim
Modified: pkg/PortfolioAnalytics/R/constrained_objective.R
===================================================================
--- pkg/PortfolioAnalytics/R/constrained_objective.R 2013-07-11 02:33:59 UTC (rev 2537)
+++ pkg/PortfolioAnalytics/R/constrained_objective.R 2013-07-11 03:40:14 UTC (rev 2538)
@@ -334,4 +334,329 @@
} else {
return(list(out=as.numeric(out),weights=w,objective_measures=tmp_return))
}
-}
\ No newline at end of file
+}
+
+#' constrained_objective_v2 2 function to calculate a numeric return value for a portfolio based on a set of constraints and objectives
+#'
+#' function to calculate a numeric return value for a portfolio based on a set of constraints,
+#' we'll try to make as few assumptions as possible, and only run objectives that are required by the user
+#'
+#' If the user has passed in either min_sum or max_sum constraints for the portfolio, or both,
+#' and are using a numerical optimization method like DEoptim, and normalize=TRUE, the default,
+#' we'll normalize the weights passed in to whichever boundary condition has been violated.
+#' If using random portfolios, all the portfolios generated will meet the constraints by construction.
+#' 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 apply the same normalization in \code{\link{optimize.portfolio}} so that the weights you see have been
+#' normalized to min_sum if the generated portfolio is smaller than min_sum or max_sum if the
+#' generated portfolio is larger than max_sum.
+#' This normalization increases the speed of optimization and convergence by several orders of magnitude in many cases.
+#'
+#' You may find that for some portfolios, normalization is not desirable, if the algorithm
+#' cannot find a direction in which to move to head towards an optimal portfolio. In these cases,
+#' it may be best to set normalize=FALSE, and penalize the portfolios if the sum of the weighting
+#' vector lies outside the min_sum and/or max_sum.
+#'
+#' Whether or not we normalize the weights using min_sum and max_sum, and are using a numerical optimization
+#' engine like DEoptim, we will penalize portfolios that violate weight constraints in much the same way
+#' we penalize other constraints. If a min_sum/max_sum normalization has not occurred, convergence
+#' can take a very long time. We currently do not allow for a non-normalized full investment constraint.
+#' Future version of this function could include this additional constraint penalty.
+#'
+#' When you are optimizing a return objective, you must specify a negative multiplier
+#' for the return objective so that the function will maximize return. If you specify a target return,
+#' any return less than your target will be penalized. If you do not specify a target return,
+#' you may need to specify a negative VTR (value to reach) , or the function will not converge.
+#' Try the maximum expected return times the multiplier (e.g. -1 or -10).
+#' Adding a return objective defaults the multiplier to -1.
+#'
+#' Additional parameters for random portfolios or \code{\link[DEoptim]{DEoptim.control}} may be passed in via \dots
+#'
+#'
+#' @param R an xts, vector, matrix, data frame, timeSeries or zoo object of asset returns
+#' @param w a vector of weights to test
+#' @param portfolio an object of type "portfolio" specifying the constraints and objectives for the optimization, see \code{\link{constraint}}
+#' @param \dots any other passthru parameters
+#' @param trace TRUE/FALSE whether to include debugging and additional detail in the output list
+#' @param normalize TRUE/FALSE whether to normalize results to min/max sum (TRUE), or let the optimizer penalize portfolios that do not conform (FALSE)
+#' @param storage TRUE/FALSE default TRUE for DEoptim with trace, otherwise FALSE. not typically user-called
+#' @seealso \code{\link{constraint}}, \code{\link{objective}}, \code{\link[DEoptim]{DEoptim.control}}
+#' @author Kris Boudt, Peter Carl, Brian G. Peterson, Ross Bennett
+#' @export
+constrained_objective_v2 <- function(w, R, portfolio, ..., trace=FALSE, normalize=TRUE, storage=FALSE)
+{
+ if (ncol(R) > length(w)) {
+ R <- R[ ,1:length(w)]
+ }
+ if(!hasArg(penalty)) penalty <- 1e4
+ N <- length(w)
+ T <- nrow(R)
+ if(hasArg(optimize_method))
+ optimize_method <- match.call(expand.dots=TRUE)$optimize_method else optimize_method <- ''
+ if(hasArg(verbose))
+ verbose <- match.call(expand.dots=TRUE)$verbose
+ else verbose <- FALSE
+
+ # get the constraints from the portfolio object
+ constraints <- get_constraints(portfolio)
+
+ # check for valid portfolio
+ if (!is.portfolio(portfolio)) {
+ stop("portfolio object passed in is not of class portfolio")
+ }
+
+ # check that the assets and the weighting vector have the same length
+ if (N != length(portfolio$assets)){
+ warning("length of portfolio asset list and weights vector do not match, results may be bogus")
+ }
+
+ out <- 0
+
+ # do the get here
+ store_output <- try(get('.objectivestorage',pos='.GlobalEnv'), silent=TRUE)
+ if(inherits(store_output,"try-error")) storage <- FALSE else storage <- TRUE
+
+ # may be replaced by fn_map later
+ if(isTRUE(normalize)){
+ 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
+ #' the normalization process produces much faster convergence,
+ #' and then we penalize parameters outside the constraints in the next block
+ if(!is.null(constraints$max_sum) & constraints$max_sum != Inf ) {
+ max_sum <- constraints$max_sum
+ if(sum(w) > max_sum) { w <- (max_sum / sum(w)) * w } # normalize to max_sum
+ }
+
+ if(!is.null(constraints$min_sum) & constraints$min_sum != -Inf ) {
+ min_sum <- constraints$min_sum
+ if(sum(w) < min_sum) { w <- (min_sum / sum(w)) * w } # normalize to min_sum
+ }
+
+ } # end min_sum and max_sum normalization
+ } else {
+ # the user wants the optimization algorithm to figure it out
+ if(!is.null(constraints$max_sum) & constraints$max_sum != Inf ) {
+ max_sum <- constraints$max_sum
+ if(sum(w) > max_sum) { out <- out + penalty * (sum(w) - max_sum) } # penalize difference to max_sum
+ }
+ if(!is.null(constraints$min_sum) & constraints$min_sum != -Inf ) {
+ min_sum <- constraints$min_sum
+ if(sum(w) < min_sum) { out <- out + penalty * (min_sum - sum(w)) } # penalize difference to min_sum
+ }
+ }
+
+ #' penalize weights outside my constraints (can be caused by normalization)
+ if (!is.null(constraints$max)){
+ max <- constraints$max
+ out <- out + sum(w[which(w > max[1:N])] - constraints$max[which(w > max[1:N])]) * penalty
+ }
+ if (!is.null(constraints$min)){
+ min <- constraints$min
+ out <- out + sum(constraints$min[which(w < min[1:N])] - w[which(w < min[1:N])]) * penalty
+ }
+
+ nargs <- list(...)
+ if(length(nargs)==0) nargs <- NULL
+ if (length('...')==0 | is.null('...')) {
+ # rm('...')
+ nargs <- NULL
+ }
+
+ nargs <- set.portfolio.moments_v2(R, portfolio, momentargs=nargs)
+
+ if(is.null(portfolio$objectives)) {
+ warning("no objectives specified in portfolio")
+ } else{
+ if(isTRUE(trace) | isTRUE(storage)) tmp_return <- list()
+ for (objective in portfolio$objectives){
+ #check for clean bits to pass in
+ if(objective$enabled){
+ tmp_measure <- NULL
+ multiplier <- objective$multiplier
+ #if(is.null(objective$arguments) | !is.list(objective$arguments)) objective$arguments<-list()
+ switch(objective$name,
+ mean =,
+ median = {
+ fun = match.fun(objective$name)
+ nargs$x <- ( R %*% w ) #do the multivariate mean/median with Kroneker product
+ },
+ sd =,
+ StdDev = {
+ fun = match.fun(StdDev)
+ },
+ mVaR =,
+ VaR = {
+ fun = match.fun(VaR)
+ if(!inherits(objective,"risk_budget_objective") & is.null(objective$arguments$portfolio_method) & is.null(nargs$portfolio_method)) nargs$portfolio_method='single'
+ if(is.null(objective$arguments$invert)) objective$arguments$invert = FALSE
+ },
+ es =,
+ mES =,
+ CVaR =,
+ cVaR =,
+ ES = {
+ fun = match.fun(ES)
+ if(!inherits(objective,"risk_budget_objective") & is.null(objective$arguments$portfolio_method)& is.null(nargs$portfolio_method)) nargs$portfolio_method='single'
+ if(is.null(objective$arguments$invert)) objective$arguments$invert = FALSE
+ },
+ turnover = {
+ fun = match.fun(turnover) # turnover function included in objectiveFUN.R
+ },
+{ # see 'S Programming p. 67 for this matching
+ fun <- try(match.fun(objective$name))
+}
+ )
+ if(is.function(fun)){
+ .formals <- formals(fun)
+ onames <- names(.formals)
+ if(is.list(objective$arguments)){
+ #TODO FIXME only do this if R and weights are in the argument list of the fn
+ if(is.null(nargs$R) | !length(nargs$R)==length(R)) nargs$R <- R
+
+ if(is.null(nargs$weights)) nargs$weights <- w
+
+ pm <- pmatch(names(objective$arguments), onames, nomatch = 0L)
+ if (any(pm == 0L))
+ warning(paste("some arguments stored for", objective$name, "do not match"))
+ # this line overwrites the names of things stored in $arguments with names from formals.
+ # I'm not sure it's a good idea, so commenting for now, until we prove we need it
+ #names(objective$arguments[pm > 0L]) <- onames[pm]
+ .formals[pm] <- objective$arguments[pm > 0L]
+ #now add dots
+ if (length(nargs)) {
+ dargs <- nargs
+ pm <- pmatch(names(dargs), onames, nomatch = 0L)
+ names(dargs[pm > 0L]) <- onames[pm]
+ .formals[pm] <- dargs[pm > 0L]
+ }
+ .formals$... <- NULL
+ }
+ } # TODO do some funky return magic here on try-error
+
+ tmp_measure <- try((do.call(fun,.formals)), silent=TRUE)
+
+ if(isTRUE(trace) | isTRUE(storage)) {
+ if(is.null(names(tmp_measure))) names(tmp_measure) <- objective$name
+ tmp_return[[objective$name]] <- tmp_measure
+ }
+
+ if(inherits(tmp_measure, "try-error")) {
+ message(paste("objective name", objective$name, "generated an error or warning:", tmp_measure))
+ next()
+ }
+
+ # now set the new value of the objective output
+ if(inherits(objective, "return_objective")){
+ if (!is.null(objective$target) & is.numeric(objective$target)){ # we have a target
+ out <- out + penalty*abs(objective$multiplier)*abs(tmp_measure - objective$target)
+ }
+ # target is null or doesn't exist, just maximize, or minimize violation of constraint
+ out <- out + objective$multiplier*tmp_measure
+ } # end handling for return objectives
+
+ if(inherits(objective, "portfolio_risk_objective")){
+ if (!is.null(objective$target) & is.numeric(objective$target)){ # we have a target
+ out <- out + penalty*abs(objective$multiplier)*abs(tmp_measure - objective$target)
+ #should we also penalize risk too low for risk targets? or is a range another objective?
+ # # half penalty for risk lower than target
+ # if( prw < (.9*Riskupper) ){ out = out + .5*(penalty*( prw - Riskupper)) }
+ }
+ # target is null or doesn't exist, just maximize, or minimize violation of constraint
+ out <- out + abs(objective$multiplier)*tmp_measure
+ } # univariate risk objectives
+
+ if(inherits(objective, "turnover_objective")){
+ if (!is.null(objective$target) & is.numeric(objective$target)){ # we have a target
+ out <- out + penalty*abs(objective$multiplier)*abs(tmp_measure - objective$target)
+ }
+ # target is null or doesn't exist, just maximize, or minimize violation of constraint
+ out <- out + abs(objective$multiplier)*tmp_measure
+ } # univariate turnover objectives
+
+ if(inherits(objective, "minmax_objective")){
+ if (!is.null(objective$min) & !is.null(objective$max)){ # we have a min and max
+ if(tmp_measure > objective$max){
+ out <- out + penalty * objective$multiplier * (tmp_measure - objective$max)
+ }
+ if(tmp_measure < objective$min){
+ out <- out + penalty * objective$multiplier * (objective$min - tmp_measure)
+ }
+ }
+ } # temporary minmax objective
+
+ if(inherits(objective, "risk_budget_objective")){
+ # setup
+
+ # out = out + penalty*sum( (percrisk-RBupper)*( percrisk > RBupper ),na.rm=TRUE ) + penalty*sum( (RBlower-percrisk)*( percrisk < RBlower ),na.rm=TRUE )
+ # add risk budget constraint
+ if(!is.null(objective$target) & is.numeric(objective$target)){
+ #in addition to a risk budget constraint, we have a univariate target
+ # the first element of the returned list is the univariate measure
+ # we'll use the univariate measure exactly like we would as a separate objective
+ out = out + penalty*abs(objective$multiplier)*abs(tmp_measure[[1]]-objective$target)
+ #should we also penalize risk too low for risk targets? or is a range another objective?
+ # # half penalty for risk lower than target
+ # if( prw < (.9*Riskupper) ){ out = out + .5*(penalty*( prw - Riskupper)) }
+ }
+ percrisk = tmp_measure[[3]] # third element is percent component contribution
+ RBupper = objective$max_prisk
+ RBlower = objective$min_prisk
+ if(!is.null(RBupper) | !is.null(RBlower)){
+ out = out + penalty * objective$multiplier * sum( (percrisk-RBupper)*( percrisk > RBupper ),na.rm=TRUE ) + penalty*sum( (RBlower-percrisk)*( percrisk < RBlower ),na.rm=TRUE )
+ }
+ # if(!is.null(objective$min_concentration)){
+ # if(isTRUE(objective$min_concentration)){
+ # max_conc<-max(tmp_measure[[2]]) #second element is the contribution in absolute terms
+ # # out=out + penalty * objective$multiplier * max_conc
+ # out = out + objective$multiplier * max_conc
+ # }
+ # }
+ # Combined min_con and min_dif to take advantage of a better concentration obj measure
+ if(!is.null(objective$min_difference) || !is.null(objective$min_concentration)){
+ if(isTRUE(objective$min_difference)){
+ # max_diff<-max(tmp_measure[[2]]-(sum(tmp_measure[[2]])/length(tmp_measure[[2]]))) #second element is the contribution in absolute terms
+ # Uses Herfindahl index to calculate concentration; added scaling perc diffs back to univariate numbers
+ max_diff <- sqrt(sum(tmp_measure[[3]]^2))/100 #third element is the contribution in percentage terms
+ # out = out + penalty * objective$multiplier * max_diff
+ out = out + penalty*objective$multiplier * max_diff
+ }
+ }
+ } # end handling of risk_budget objective
+
+ } # end enabled check
+ } # end loop over objectives
+ } # end objectives processing
+
+ if(isTRUE(verbose)) {
+ print('weights: ')
+ print(paste(w,' '))
+ print(paste("output of objective function", out))
+ print(unlist(tmp_return))
+ }
+
+ if(is.na(out) | is.nan(out) | is.null(out)){
+ #this should never happen
+ warning('NA or NaN produced in objective function for weights ',w)
+ out <- penalty
+ }
+
+ #return
+ if (isTRUE(storage)){
+ #add the new objective results
+ store_output[[length(store_output)+1]] <- list(out=as.numeric(out), weights=w, objective_measures=tmp_return)
+ # do the assign here
+ assign('.objectivestorage', store_output, pos='.GlobalEnv')
+ }
+ if(!isTRUE(trace)){
+ return(out)
+ } else {
+ return(list(out=as.numeric(out), weights=w, objective_measures=tmp_return))
+ }
+}
Modified: pkg/PortfolioAnalytics/R/optimize.portfolio.R
===================================================================
--- pkg/PortfolioAnalytics/R/optimize.portfolio.R 2013-07-11 02:33:59 UTC (rev 2537)
+++ pkg/PortfolioAnalytics/R/optimize.portfolio.R 2013-07-11 03:40:14 UTC (rev 2538)
@@ -476,9 +476,9 @@
momentFUN='set.portfolio.moments_v2'
)
{
- optimize_method=optimize_method[1]
- tmptrace=NULL
- start_t<-Sys.time()
+ optimize_method <- optimize_method[1]
+ tmptrace <- NULL
+ start_t <- Sys.time()
#store the call for later
call <- match.call()
@@ -500,6 +500,9 @@
dotargs <- list(...)
+ # Get the constraints from the portfolio object
+ constraints <- get_constraints(portfolio)
+
# set portfolio moments only once
if(!is.function(momentFUN)){
momentFUN <- match.fun(momentFUN)
@@ -515,6 +518,131 @@
} else {
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(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 else rpseed=TRUE
+ if(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
+ # This part should still work, but will need to change random_portfolios over to accept portfolio object
+ 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(rpconstraints=rpconstraint, permutations=NP)
+ DEcformals$initialpop=rp
+ }
+ controlDE <- do.call(DEoptim.control, DEcformals)
+
+ # need to modify constrained_objective to accept a portfolio object
+ minw = try(DEoptim( constrained_objective_v2, lower=lower[1:N], upper=upper[1:N], control=controlDE, R=R, portfolio=portfolio, 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_v2(w=weights, R=R, portfolio, 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
+
+ # 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))
+ message(c("elapsed time:", end_t-start_t))
+ out$portfolio <- portfolio
+ 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)
}
#' portfolio optimization with support for rebalancing or rolling periods
More information about the Returnanalytics-commits
mailing list