[Returnanalytics-commits] r3299 - in pkg/PortfolioAnalytics: . R demo

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Sat Jan 11 22:23:39 CET 2014


Author: rossbennett34
Date: 2014-01-11 22:23:39 +0100 (Sat, 11 Jan 2014)
New Revision: 3299

Added:
   pkg/PortfolioAnalytics/demo/demo_roi_solvers.R
Modified:
   pkg/PortfolioAnalytics/DESCRIPTION
   pkg/PortfolioAnalytics/R/optFUN.R
   pkg/PortfolioAnalytics/R/optimize.portfolio.R
   pkg/PortfolioAnalytics/demo/00Index
Log:
Minor modifications to optFUN and optimize.portfolio for easy way to specify solver for ROI per discussion with Doug.

Modified: pkg/PortfolioAnalytics/DESCRIPTION
===================================================================
--- pkg/PortfolioAnalytics/DESCRIPTION	2014-01-10 23:43:31 UTC (rev 3298)
+++ pkg/PortfolioAnalytics/DESCRIPTION	2014-01-11 21:23:39 UTC (rev 3299)
@@ -24,6 +24,7 @@
     ROI (>= 0.1.0),
     ROI.plugin.glpk (>= 0.0.2),
     ROI.plugin.quadprog (>= 0.0.2),
+    ROI.plugin.symphony (>= 0.0.2),
     pso,
     GenSA,
     corpcor

Modified: pkg/PortfolioAnalytics/R/optFUN.R
===================================================================
--- pkg/PortfolioAnalytics/R/optFUN.R	2014-01-10 23:43:31 UTC (rev 3298)
+++ pkg/PortfolioAnalytics/R/optFUN.R	2014-01-11 21:23:39 UTC (rev 3299)
@@ -1,8 +1,9 @@
 
 ##### GMV and QU QP Function #####
-#' Optimization function to solve minimum variance or maximum quadratic utility problems
+#' GMV/QU QP Optimization
 #' 
-#' This function is called by optimize.portfolio to solve minimum variance or maximum quadratic utility problems
+#' This function is called by optimize.portfolio to solve minimum variance or 
+#' maximum quadratic utility problems
 #' 
 #' @param R xts object of asset returns
 #' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
@@ -10,9 +11,10 @@
 #' @param lambda risk_aversion parameter
 #' @param target target return value
 #' @param lambda_hhi concentration aversion parameter
-#' @param conc_groups list of vectors specifying the groups of the assets. 
+#' @param conc_groups list of vectors specifying the groups of the assets.
+#' @param solver solver to use
 #' @author Ross Bennett
-gmv_opt <- function(R, constraints, moments, lambda, target, lambda_hhi, conc_groups){
+gmv_opt <- function(R, constraints, moments, lambda, target, lambda_hhi, conc_groups, solver="quadprog"){
   stopifnot("package:ROI" %in% search() || require("ROI", quietly = TRUE))
   stopifnot("package:ROI.plugin.quadprog" %in% search() || require("ROI.plugin.quadprog", quietly = TRUE))
   
@@ -119,7 +121,7 @@
   # set up the optimization problem and solve
   opt.prob <- OP(objective=ROI_objective, 
                        constraints=L_constraint(L=Amat, dir=dir.vec, rhs=rhs.vec))
-  result <- ROI_solve(x=opt.prob, solver="quadprog")
+  result <- ROI_solve(x=opt.prob, solver=solver)
   
   # result <- try(solve.QP(Dmat=Dmat, dvec=dvec, Amat=t(Amat), bvec=rhs.vec, meq=meq), silent=TRUE)
   if(inherits(x=result, "try-error")) stop(paste("No solution found:", result))
@@ -152,16 +154,17 @@
 }
 
 ##### Maximize Return LP Function #####
-#' Optimization function to solve minimum variance or maximum quadratic utility problems
+#' Maximum Return LP Optimization
 #' 
-#' This function is called by optimize.portfolio to solve minimum variance or maximum quadratic utility problems
+#' This function is called by optimize.portfolio to solve maximum return
 #' 
 #' @param R xts object of asset returns
 #' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
 #' @param moments object of moments computed based on objective functions
 #' @param target target return value
+#' @param solver solver to use
 #' @author Ross Bennett
-maxret_opt <- function(R, moments, constraints, target){
+maxret_opt <- function(R, moments, constraints, target, solver="glpk"){
   stopifnot("package:ROI" %in% search() || require("ROI",quietly = TRUE))
   stopifnot("package:ROI.plugin.glpk" %in% search() || require("ROI.plugin.glpk",quietly = TRUE))
   
@@ -224,7 +227,7 @@
   opt.prob <- OP(objective=ROI_objective, 
                  constraints=L_constraint(L=Amat, dir=dir.vec, rhs=rhs.vec),
                  bounds=bnds)
-  roi.result <- ROI_solve(x=opt.prob, solver="glpk")
+  roi.result <- ROI_solve(x=opt.prob, solver=solver)
   
   # roi.result <- Rglpk_solve_LP(obj=objL, mat=Amat, dir=dir.vec, rhs=rhs.vec, bounds=bnds)
   
@@ -255,16 +258,18 @@
 }
 
 ##### Maximize Return MILP Function #####
-#' Optimization function to solve maximum return problems
+#' Maximum Return MILP Optimization
 #' 
-#' This function is called by optimize.portfolio to solve maximum return problems via mixed integer linear programming.
+#' This function is called by optimize.portfolio to solve maximum return 
+#' problems via mixed integer linear programming.
 #' 
 #' @param R xts object of asset returns
 #' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
 #' @param moments object of moments computed based on objective functions
 #' @param target target return value
+#' @param solver solver to use
 #' @author Ross Bennett
-maxret_milp_opt <- function(R, constraints, moments, target){
+maxret_milp_opt <- function(R, constraints, moments, target, solver="glpk"){
   stopifnot("package:ROI" %in% search() || require("ROI",quietly = TRUE))
   stopifnot("package:ROI.plugin.glpk" %in% search() || require("ROI.plugin.glpk",quietly = TRUE))
   
@@ -356,7 +361,7 @@
   opt.prob <- OP(objective=ROI_objective, 
                  constraints=L_constraint(L=Amat, dir=dir, rhs=rhs),
                  bounds=bnds, types=types)
-  roi.result <- try(ROI_solve(x=opt.prob, solver="glpk"), silent=TRUE)
+  roi.result <- try(ROI_solve(x=opt.prob, solver=solver), silent=TRUE)
   if(inherits(roi.result, "try-error")) stop(paste("No solution found:", roi.result))
   
   # Weights
@@ -377,7 +382,7 @@
 }
 
 ##### Minimize ETL LP Function #####
-#' Optimization function to solve minimum ETL problems
+#' Minimum ETL LP Optimization
 #' 
 #' This function is called by optimize.portfolio to solve minimum ETL problems.
 #' 
@@ -386,8 +391,9 @@
 #' @param moments object of moments computed based on objective functions
 #' @param target target return value
 #' @param alpha alpha value for ETL/ES/CVaR
+#' @param solver solver to use
 #' @author Ross Bennett
-etl_opt <- function(R, constraints, moments, target, alpha){
+etl_opt <- function(R, constraints, moments, target, alpha, solver="glpk"){
   stopifnot("package:ROI" %in% search() || require("ROI",quietly = TRUE))
   stopifnot("package:ROI.plugin.glpk" %in% search() || require("ROI.plugin.glpk",quietly = TRUE))
   
@@ -442,7 +448,7 @@
   opt.prob <- OP(objective=ROI_objective, 
                        constraints=L_constraint(L=Amat, dir=dir.vec, rhs=rhs.vec),
                        bounds=bnds)
-  roi.result <- try(ROI_solve(x=opt.prob, solver="glpk"), silent=TRUE)
+  roi.result <- try(ROI_solve(x=opt.prob, solver=solver), silent=TRUE)
   if(inherits(x=roi.result, "try-error")) stop(paste("No solution found:", roi.result))
   
   weights <- roi.result$solution[1:N]
@@ -473,17 +479,19 @@
 }
 
 ##### Minimize ETL MILP Function #####
-#' Optimization function to solve minimum ETL problems
+#' Minimum ETL MILP Optimization
 #' 
-#' This function is called by optimize.portfolio to solve minimum ETL problems via mixed integer linear programming.
+#' This function is called by optimize.portfolio to solve minimum ETL problems 
+#' via mixed integer linear programming.
 #' 
 #' @param R xts object of asset returns
 #' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
 #' @param moments object of moments computed based on objective functions
 #' @param target target return value
 #' @param alpha alpha value for ETL/ES/CVaR
+#' @param solver solver to use
 #' @author Ross Bennett
-etl_milp_opt <- function(R, constraints, moments, target, alpha){
+etl_milp_opt <- function(R, constraints, moments, target, alpha, solver="glpk"){
   stopifnot("package:ROI" %in% search() || require("ROI",quietly = TRUE))
   stopifnot("package:ROI.plugin.glpk" %in% search() || require("ROI.plugin.glpk",quietly = TRUE))
   
@@ -589,7 +597,7 @@
   opt.prob <- OP(objective=ROI_objective, 
                  constraints=L_constraint(L=tmpAmat, dir=dir, rhs=rhs),
                  bounds=bnds, types=types)
-  roi.result <- ROI_solve(x=opt.prob, solver="glpk")
+  roi.result <- ROI_solve(x=opt.prob, solver=solver)
   
   # The Rglpk solvers status returns an an integer with status information
   # about the solution returned: 0 if the optimal solution was found, a 
@@ -626,9 +634,10 @@
 }
 
 ##### minimize variance or maximize quadratic utility with turnover constraints #####
-#' Optimization function to solve minimum variance or maximum quadratic utility problems with turnover constraint
+#' GMV/QU QP Optimization with Turnover Constraint
 #' 
-#' This function is called by optimize.portfolio to solve minimum variance or maximum quadratic utility problems
+#' This function is called by optimize.portfolio to solve minimum variance or 
+#' maximum quadratic utility problems with turnover constraint
 #' 
 #' @param R xts object of asset returns
 #' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
@@ -636,8 +645,9 @@
 #' @param lambda risk_aversion parameter
 #' @param target target return value
 #' @param init_weights initial weights to compute turnover
+#' @param solver solver to use
 #' @author Ross Bennett
-gmv_opt_toc <- function(R, constraints, moments, lambda, target, init_weights){
+gmv_opt_toc <- function(R, constraints, moments, lambda, target, init_weights, solver="quadprog"){
   # function for minimum variance or max quadratic utility problems
   stopifnot("package:corpcor" %in% search() || require("corpcor",quietly = TRUE))
   stopifnot("package:ROI" %in% search() || require("ROI", quietly = TRUE))
@@ -748,7 +758,7 @@
   opt.prob <- OP(objective=ROI_objective, 
                  constraints=L_constraint(L=Amat, dir=dir, rhs=rhs))
   
-  roi.result <- try(ROI_solve(x=opt.prob, solver="quadprog"), silent=TRUE)
+  roi.result <- try(ROI_solve(x=opt.prob, solver=solver), silent=TRUE)
   
   if(inherits(roi.result, "try-error")) stop(paste("No solution found:", roi.result))
   
@@ -779,8 +789,21 @@
   return(out)
 }
 
-# proportional transaction cost constraint
-gmv_opt_ptc <- function(R, constraints, moments, lambda, target, init_weights){
+##### minimize variance or maximize quadratic utility with proportional transactioncosts constraints #####
+#' GMV/QU QP Optimization with Proportional Transaction Cost Constraint
+#' 
+#' This function is called by optimize.portfolio to solve minimum variance or 
+#' maximum quadratic utility problems with proportional transaction cost constraint
+#' 
+#' @param R xts object of asset returns
+#' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
+#' @param moments object of moments computed based on objective functions
+#' @param lambda risk_aversion parameter
+#' @param target target return value
+#' @param init_weights initial weights to compute turnover
+#' @param solver solver to use
+#' @author Ross Bennett
+gmv_opt_ptc <- function(R, constraints, moments, lambda, target, init_weights, solver="quadprog"){
   # function for minimum variance or max quadratic utility problems
   # modifying ProportionalCostOpt function from MPO package
   stopifnot("package:corpcor" %in% search() || require("corpcor", quietly = TRUE))
@@ -880,7 +903,7 @@
   
   opt.prob <- OP(objective=ROI_objective, 
                  constraints=L_constraint(L=Amat, dir=dir, rhs=rhs))
-  roi.result <- try(ROI_solve(x=opt.prob, solver="quadprog"), silent=TRUE)
+  roi.result <- try(ROI_solve(x=opt.prob, solver=solver), silent=TRUE)
 
   if(inherits(roi.result, "try-error")) stop(paste("No solution found:", roi.result))
   
@@ -915,7 +938,7 @@
 }
 
 
-mean_etl_opt <- function(R, constraints, moments, target, alpha, tol=.Machine$double.eps^0.5, maxit=50){
+mean_etl_opt <- function(R, constraints, moments, target, alpha, solver="glpk", tol=.Machine$double.eps^0.5, maxit=50){
   # This function returns the target mean return that maximizes mean / etl (i.e. starr)
   
   # if all(moments$mean == 0) then the user did not specify mean as an objective,
@@ -929,17 +952,17 @@
   
   # Find the maximum return
   if(!is.null(constraints$max_pos)){
-    max_ret <- maxret_milp_opt(R=R, constraints=constraints, moments=moments, target=NA)
+    max_ret <- maxret_milp_opt(R=R, constraints=constraints, moments=moments, target=NA, solver=solver)
   } else {
-    max_ret <- maxret_opt(R=R, moments=moments, constraints=constraints, target=NA)
+    max_ret <- maxret_opt(R=R, moments=moments, constraints=constraints, target=NA, solver=solver)
   }
   max_mean <- as.numeric(-max_ret$out)
   
   # Find the starr at the maximum etl portfolio
   if(!is.null(constraints$max_pos)){
-    ub_etl <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=max_mean, alpha=alpha)
+    ub_etl <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=max_mean, alpha=alpha, solver=solver)
   } else {
-    ub_etl <- etl_opt(R=R, constraints=constraints, moments=moments, target=max_mean, alpha=alpha)
+    ub_etl <- etl_opt(R=R, constraints=constraints, moments=moments, target=max_mean, alpha=alpha, solver=solver)
   }
   ub_weights <- matrix(ub_etl$weights, ncol=1)
   ub_mean <- as.numeric(t(ub_weights) %*% fmean)
@@ -954,9 +977,9 @@
   
   # Find the starr at the minimum etl portfolio
   if(!is.null(constraints$max_pos)){
-    lb_etl <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=NA, alpha=alpha)
+    lb_etl <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=NA, alpha=alpha, solver=solver)
   } else {
-    lb_etl <- etl_opt(R=R, constraints=constraints, moments=moments, target=NA, alpha=alpha)
+    lb_etl <- etl_opt(R=R, constraints=constraints, moments=moments, target=NA, alpha=alpha, solver=solver)
   }
   lb_weights <- matrix(lb_etl$weights)  
   lb_mean <- as.numeric(t(lb_weights) %*% fmean)  
@@ -986,9 +1009,9 @@
     # Find the starr at the mean return midpoint
     new_ret <- (lb_mean + ub_mean) / 2
     if(!is.null(constraints$max_pos)){
-      mid <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha)
+      mid <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha, solver=solver)
     } else {
-      mid <- etl_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha)
+      mid <- etl_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha, solver=solver)
     }
     # print(mid)
     mid_weights <- matrix(mid$weights, ncol=1)
@@ -1009,9 +1032,9 @@
       ub_starr <- mid_starr
       new_ret <- (lb_mean + ub_mean) / 2
       if(!is.null(constraints$max_pos)){
-        mid <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha)
+        mid <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha, solver=solver)
       } else {
-        mid <- etl_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha)
+        mid <- etl_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha, solver=solver)
       }
       mid_weights <- matrix(mid$weights, ncol=1)
       mid_mean <- as.numeric(t(mid_weights) %*% fmean)
@@ -1025,9 +1048,9 @@
       lb_starr <- mid_starr
       new_ret <- (lb_mean + ub_mean) / 2
       if(!is.null(constraints$max_pos)){
-        mid <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha)
+        mid <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha, solver=solver)
       } else {
-        mid <- etl_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha)
+        mid <- etl_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha, solver=solver)
       }
       mid_weights <- matrix(mid$weights, ncol=1)
       mid_mean <- as.numeric(t(mid_weights) %*% fmean)
@@ -1041,7 +1064,7 @@
   return(new_ret)
 }
 
-max_sr_opt <- function(R, constraints, moments, lambda, target, lambda_hhi, conc_groups, tol=.Machine$double.eps^0.5, maxit=50){
+max_sr_opt <- function(R, constraints, moments, lambda, target, lambda_hhi, conc_groups, solver="quadprog", tol=.Machine$double.eps^0.5, maxit=50){
   # This function returns the target mean return that maximizes mean / sd (i.e. sharpe ratio)
   
   # get the forecast mean from moments
@@ -1061,7 +1084,7 @@
   # Calculate the sr at the miminum var portfolio
   tmpmoments <- moments
   tmpmoments$mean <- rep(0, length(moments$mean))
-  lb_sr <- gmv_opt(R=R, constraints=constraints, moments=tmpmoments, lambda=1, target=NA, lambda_hhi=lambda_hhi, conc_groups=conc_groups)
+  lb_sr <- gmv_opt(R=R, constraints=constraints, moments=tmpmoments, lambda=1, target=NA, lambda_hhi=lambda_hhi, conc_groups=conc_groups, solver=solver)
   lb_weights <- matrix(lb_sr$weights)
   lb_mean <- as.numeric(t(lb_weights) %*% fmean)
   lb_sd <- as.numeric(sqrt(t(lb_weights) %*% moments$var %*% lb_weights))
@@ -1079,7 +1102,7 @@
     
     # Find the starr at the mean return midpoint
     new_ret <- (lb_mean + ub_mean) / 2
-    mid <- gmv_opt(R=R, constraints=constraints, moments=tmpmoments, lambda=1, target=new_ret, lambda_hhi=lambda_hhi, conc_groups=conc_groups)
+    mid <- gmv_opt(R=R, constraints=constraints, moments=tmpmoments, lambda=1, target=new_ret, lambda_hhi=lambda_hhi, conc_groups=conc_groups, solver=solver)
     mid_weights <- matrix(mid$weights, ncol=1)
     mid_mean <- as.numeric(t(mid_weights) %*% fmean)
     mid_sd <- as.numeric(sqrt(t(mid_weights) %*% moments$var %*% mid_weights))
@@ -1096,7 +1119,7 @@
       ub_mean <- mid_mean
       ub_sr <- mid_sr
       new_ret <- (lb_mean + ub_mean) / 2
-      mid <- gmv_opt(R=R, constraints=constraints, moments=tmpmoments, lambda=1, target=new_ret, lambda_hhi=lambda_hhi, conc_groups=conc_groups)
+      mid <- gmv_opt(R=R, constraints=constraints, moments=tmpmoments, lambda=1, target=new_ret, lambda_hhi=lambda_hhi, conc_groups=conc_groups, solver=solver)
       mid_weights <- matrix(mid$weights, ncol=1)
       mid_mean <- as.numeric(t(mid_weights) %*% fmean)
       mid_sd <- as.numeric(sqrt(t(mid_weights) %*% moments$var %*% mid_weights))
@@ -1106,7 +1129,7 @@
       lb_mean <- mid_mean
       lb_sr <- mid_sr
       new_ret <- (lb_mean + ub_mean) / 2
-      mid <- gmv_opt(R=R, constraints=constraints, moments=tmpmoments, lambda=1, target=new_ret, lambda_hhi=lambda_hhi, conc_groups=conc_groups)
+      mid <- gmv_opt(R=R, constraints=constraints, moments=tmpmoments, lambda=1, target=new_ret, lambda_hhi=lambda_hhi, conc_groups=conc_groups, solver=solver)
       mid_weights <- matrix(mid$weights, ncol=1)
       mid_mean <- as.numeric(t(mid_weights) %*% fmean)
       mid_sd <- as.numeric(sqrt(t(mid_weights) %*% moments$var %*% mid_weights))

Modified: pkg/PortfolioAnalytics/R/optimize.portfolio.R
===================================================================
--- pkg/PortfolioAnalytics/R/optimize.portfolio.R	2014-01-10 23:43:31 UTC (rev 3298)
+++ pkg/PortfolioAnalytics/R/optimize.portfolio.R	2014-01-11 21:23:39 UTC (rev 3299)
@@ -438,7 +438,7 @@
   portfolio=NULL,
   constraints=NULL,
   objectives=NULL,
-  optimize_method=c("DEoptim","random","ROI","ROI_old","pso","GenSA"),
+  optimize_method=c("DEoptim","random","ROI","pso","GenSA"),
   search_size=20000,
   trace=FALSE, ...,
   rp=NULL,
@@ -725,11 +725,11 @@
     
   } ## end case for random
   
-  if(optimize_method == "ROI"){
-    # This takes in a regular portfolio object and extracts the desired business objectives
-    # and converts them to matrix form to be inputed into a closed form solver
-    # retrieve the objectives to minimize, these should either be "var" and/or "mean"
-    # we can either miniminze variance or maximize quiadratic utility (we will be minimizing the neg. quad. utility)
+  roi_solvers <- c("ROI", "quadprog", "glpk", "symphony", "ipop", "cplex")
+  if(optimize_method %in% roi_solvers){
+    # 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)){
@@ -737,13 +737,6 @@
     } else {
       target <- NA
     }
-    # comment out so concentration aversion can only be specified as an objective
-    # because it is added to the quadratic objective term for QP problems (minvar and qu)
-    # if(!is.null(constraints$conc_aversion)){
-    #  lambda_hhi <- constraints$conc_aversion
-    #} else {
-    #  lambda_hhi <- 0
-    #}
     lambda <- 1
     
     # list of valid objective names for ROI solvers
@@ -802,6 +795,13 @@
     }
     
     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)){
@@ -810,14 +810,14 @@
           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)
+          qp_result <- gmv_opt_toc(R=R, constraints=constraints, moments=moments, lambda=lambda, target=target, init_weights=portfolio$assets, solver=solver)
           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)
+          qp_result <- gmv_opt_ptc(R=R, constraints=constraints, moments=moments, lambda=lambda, target=target, init_weights=portfolio$assets, solver=solver)
           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
@@ -827,12 +827,12 @@
         # 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=lambda, target=target, lambda_hhi=lambda_hhi, conc_groups=conc_groups)
+          target <- max_sr_opt(R=R, constraints=constraints, moments=moments, lambda=lambda, target=target, lambda_hhi=lambda_hhi, conc_groups=conc_groups, solver=solver)
           # 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)
+        roi_result <- gmv_opt(R=R, constraints=constraints, moments=moments, lambda=lambda, target=target, lambda_hhi=lambda_hhi, conc_groups=conc_groups, solver=solver)
         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
@@ -846,17 +846,25 @@
       }
     }
     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)
+        roi_result <- maxret_milp_opt(R=R, constraints=constraints, moments=moments, target=target, solver=solver)
         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)
+        roi_result <- maxret_opt(R=R, constraints=constraints, moments=moments, target=target, solver=solver)
         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
@@ -864,6 +872,14 @@
       }
     }
     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
@@ -872,12 +888,12 @@
       # 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, target=target, alpha=alpha)
+        target <- mean_etl_opt(R=R, constraints=constraints, moments=moments, target=target, alpha=alpha, solver=solver)
         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)
+        roi_result <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=target, alpha=alpha, solver=solver)
         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
@@ -888,7 +904,7 @@
         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)
+        roi_result <- etl_opt(R=R, constraints=constraints, moments=moments, target=target, alpha=alpha, solver=solver)
         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
@@ -898,6 +914,8 @@
         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
@@ -1032,6 +1050,13 @@
 #' 
 #' 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.}
@@ -1063,7 +1088,8 @@
 #' @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","ROI_old", "pso", "GenSA".  For using \code{ROI_old}, need to use a constraint_ROI object in constraints. For using \code{ROI}, pass standard \code{constratint} object in \code{constraints} argument.  Presently, ROI has plugins for \code{quadprog} and \code{Rglpk}.
+#' @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

Modified: pkg/PortfolioAnalytics/demo/00Index
===================================================================
--- pkg/PortfolioAnalytics/demo/00Index	2014-01-10 23:43:31 UTC (rev 3298)
+++ pkg/PortfolioAnalytics/demo/00Index	2014-01-11 21:23:39 UTC (rev 3299)
@@ -24,4 +24,5 @@
 demo_min_StdDev Demonstrate objective to minimize portfolio standard deviation.
 demo_min_expected_shortfall Demonstrate objective to minimize expected shortfall.
 demo_risk_budgets Demonstrate using risk budget objectives.
+demo_roi_solvers Demonstrate specifying a solver using ROI
 

Added: pkg/PortfolioAnalytics/demo/demo_roi_solvers.R
===================================================================
--- pkg/PortfolioAnalytics/demo/demo_roi_solvers.R	                        (rev 0)
+++ pkg/PortfolioAnalytics/demo/demo_roi_solvers.R	2014-01-11 21:23:39 UTC (rev 3299)
@@ -0,0 +1,53 @@
+
+library(PortfolioAnalytics)
+library(ROI)
+library(ROI.plugin.quadprog)
+library(ROI.plugin.glpk)
+library(ROI.plugin.symphony)
+
+data(edhec)
+R <- edhec[, 1:4]
+funds <- colnames(R)
+
+# Set up initial portfolio with basic constraints
+init.portf <- portfolio.spec(assets=funds)
+init.portf <- add.constraint(portfolio=init.portf, type="full_investment")
+init.portf <- add.constraint(portfolio=init.portf, type="long_only")
+
+# Add objectives
+maxret.portf <- add.objective(portfolio=init.portf, type="return", name="mean")
+mines.portf <- add.objective(portfolio=init.portf, type="risk", name="ES")
+minsd.portf <- add.objective(portfolio=init.portf, type="risk", name="StdDev")
+qu.portf <- add.objective(portfolio=init.portf, type="risk", name="StdDev", 
+                          risk_aversion=0.25)
+qu.portf <- add.objective(portfolio=qu.portf, type="return", name="mean")
+
+# Maximize return
+opt.maxret.roi <- optimize.portfolio(R, maxret.portf, optimize_method="ROI")
+opt.maxret.glpk <- optimize.portfolio(R, maxret.portf, optimize_method="glpk")
+opt.maxret.symphony <- optimize.portfolio(R, maxret.portf, optimize_method="symphony")
+all.equal(extractStats(opt.maxret.roi), extractStats(opt.maxret.glpk))
+all.equal(extractStats(opt.maxret.roi), extractStats(opt.maxret.symphony))
+# This fails because an optimization problem with a linear objective cannot
+# be solved with a quadratic programming solver
+# opt.maxret.qp <- optimize.portfolio(R, maxret.portf, optimize_method="quadprog")
+
+# Minimize ES
+opt.mines.roi <- optimize.portfolio(R, mines.portf, optimize_method="ROI")
+opt.mines.glpk <- optimize.portfolio(R, mines.portf, optimize_method="glpk")
+opt.mines.symphony <- optimize.portfolio(R, mines.portf, optimize_method="symphony")
+all.equal(extractStats(opt.mines.roi), extractStats(opt.mines.glpk))
+all.equal(extractStats(opt.mines.roi), extractStats(opt.mines.symphony))
+
+# Minimize standard deviation
+opt.minsd.roi <- optimize.portfolio(R, minsd.portf, optimize_method="ROI")
+opt.minsd.qp <- optimize.portfolio(R, minsd.portf, optimize_method="quadprog")
+all.equal(extractStats(opt.minsd.roi), extractStats(opt.minsd.qp))
+# This fails because an optimization problem with a quadratic objective cannot
+# be solved with a linear programming solver
+# opt.minsd.glpk <- optimize.portfolio(R, minsd.portf, optimize_method="glpk")
+
+# Maximize quadratic utility
[TRUNCATED]

To get the complete diff run:
    svnlook diff /svnroot/returnanalytics -r 3299


More information about the Returnanalytics-commits mailing list