CoCalc -- optimize.portfolio.R

GitHub Repository: braverock/portfolioanalytics
Path: blob/master/R/optimize.portfolio.R
¹⁴³³ views
1
###############################################################################
2
# R (https://r-project.org/) Numeric Methods for Optimization of Portfolios
3
#
4
# Copyright (c) 2004-2023 Brian G. Peterson, Peter Carl, Ross Bennett, Kris Boudt, Xiaokang Feng, Xinran Zhao
5
#
6
# This library is distributed under the terms of the GNU Public License (GPL)
7
# for full details see the file COPYING
8
#
9
# $Id$
10
#
11
###############################################################################
12

13

14
#' @rdname optimize.portfolio
15
#' @name optimize.portfolio
16
#' @export optimize.portfolio_v1
17
optimize.portfolio_v1 <- function(
18
		R,
19
		constraints,
20
		optimize_method=c("DEoptim","random","ROI","ROI_old","pso","GenSA"), 
21
		search_size=20000, 
22
		trace=FALSE, ..., 
23
		rp=NULL,
24
		momentFUN='set.portfolio.moments_v1'
25
)
26
{
27
  optimize_method=optimize_method[1]
28
  tmptrace=NULL
29
  start_t<-Sys.time()
30

31
  #store the call for later
32
  call <- match.call()
33

34
  if (is.null(constraints) | !is.constraint(constraints)){
35
      stop("you must pass in an object of class constraints to control the optimization")
36
  }
37
  
38
  R <- checkData(R)
39
  N = length(constraints$assets)
40
  if (ncol(R)>N) {
41
      R=R[,names(constraints$assets)]
42
  }
43
  T = nrow(R)
44
    
45
  out=list()
46
  
47
  weights=NULL
48
    
49
  dotargs <-list(...)    
50
  
51
  # set portfolio moments only once
52
  if(!is.function(momentFUN)){
53
	  momentFUN<-match.fun(momentFUN)
54
  }	
55
  # TODO FIXME should match formals later
56
  #dotargs <- set.portfolio.moments(R, constraints, momentargs=dotargs)
57
  .mformals <- dotargs
58
  .mformals$R <- R
59
  .mformals$constraints <- constraints
60
  mout <- try((do.call(momentFUN,.mformals)) ,silent=TRUE)	
61
  if(inherits(mout,"try-error")) { 
62
	  message(paste("portfolio moment function failed with message",mout))
63
  } else {
64
	  dotargs <- c(dotargs,mout)
65
  }
66
	  
67
  normalize_weights <- function(weights){
68
      # normalize results if necessary
69
      if(!is.null(constraints$min_sum) | !is.null(constraints$max_sum)){
70
          # the user has passed in either min_sum or max_sum constraints for the portfolio, or both.
71
          # we'll normalize the weights passed in to whichever boundary condition has been violated
72
          # NOTE: this means that the weights produced by a numeric optimization algorithm like DEoptim
73
          # might violate your constraints, so you'd need to renormalize them after optimizing
74
          # we'll create functions for that so the user is less likely to mess it up.
75
          
76
          # NOTE: need to normalize in the optimization wrapper too before we return, since we've normalized in here
77
          # In Kris' original function, this was manifested as a full investment constraint
78
          if(!is.null(constraints$max_sum) & constraints$max_sum != Inf ) {
79
              max_sum=constraints$max_sum
80
              if(sum(weights)>max_sum) { weights<-(max_sum/sum(weights))*weights } # normalize to max_sum
81
          }
82
          
83
          if(!is.null(constraints$min_sum) & constraints$min_sum != -Inf ) {
84
              min_sum=constraints$min_sum
85
              if(sum(weights)<min_sum) { weights<-(min_sum/sum(weights))*weights } # normalize to min_sum
86
          }
87
          
88
      } # end min_sum and max_sum normalization
89
      return(weights)
90
  }
91
  
92
  if(optimize_method=="DEoptim"){
93
    stopifnot("package:DEoptim" %in% search()  ||  requireNamespace("DEoptim",quietly = TRUE) )
94
    # DEoptim does 200 generations by default, so lets set the size of each generation to search_size/200)
95
    if(hasArg(itermax)) itermax=match.call(expand.dots=TRUE)$itermax else itermax=N*50
96
    NP = round(search_size/itermax)
97
    if(NP<(N*10)) NP <- N*10
98
    if(NP>2000) NP=2000
99
    if(!hasArg(itermax)) {
100
        itermax<-round(search_size/NP)
101
        if(itermax<50) itermax=50 #set minimum number of generations
102
    }
103
    
104
    #check to see whether we need to disable foreach for parallel optimization, esp if called from inside foreach
105
    if(hasArg(parallel)) parallel=match.call(expand.dots=TRUE)$parallel else parallel=TRUE
106
    if(!isTRUE(parallel) && 'package:foreach' %in% search()){
107
        foreach::registerDoSEQ()
108
    }
109
    
110
    DEcformals  <- formals(DEoptim::DEoptim.control)
111
    DEcargs <- names(DEcformals)
112
    if( is.list(dotargs) ){
113
        pm <- pmatch(names(dotargs), DEcargs, nomatch = 0L)
114
        names(dotargs[pm > 0L]) <- DEcargs[pm]
115
        DEcformals$NP <- NP
116
        DEcformals$itermax <- itermax
117
        DEcformals[pm] <- dotargs[pm > 0L]
118
		if(!hasArg(strategy)) {
119
		  # use DE/current-to-p-best/1
120
		  strategy=6
121
      DEcformals$strategy=strategy
122
      }
123
		if(!hasArg(reltol)) {
124
		  # 1/1000 of 1% change in objective is significant
125
		  reltol=.000001
126
      DEcformals$reltol=reltol
127
      }
128
		if(!hasArg(steptol)) {
129
		  # number of assets times 1.5 tries to improve
130
		  steptol=round(N*1.5)
131
      DEcformals$steptol=steptol
132
      }
133
		if(!hasArg(c)) {
134
		  # JADE mutation parameter, this could maybe use some adjustment
135
		  tmp.c=.4
136
      DEcformals$c=tmp.c
137
      }
138
        if(!hasArg(storepopfrom)) {
139
          storepopfrom=1
140
          DEcformals$storepopfrom=storepopfrom
141
        }
142
        if(isTRUE(parallel) && 'package:foreach' %in% search()){
143
            if(!hasArg(parallelType)) {
144
              #use all cores
145
              parallelType=2
146
              DEcformals$parallelType=parallelType
147
              }
148
            if(!hasArg(packages)) {
149
              #use all packages
150
              packages <- names(sessionInfo()$otherPkgs)
151
              DEcformals$packages <- packages
152
              }
153
        }
154
		 
155
        #TODO FIXME also check for a passed in controlDE list, including checking its class, and match formals
156
    }
157
    
158
    if(isTRUE(trace)) { 
159
        #we can't pass trace=TRUE into constrained objective with DEoptim, because it expects a single numeric return
160
        tmptrace=trace 
161
        assign('.objectivestorage', list(), as.environment(.storage))
162
        trace=FALSE
163
    } 
164
    
165
    # get upper and lower weights parameters from constraints
166
    upper = constraints$max
167
    lower = constraints$min
168

169
    if(hasArg(rpseed)){ 
170
      seed <- match.call(expand.dots=TRUE)$rpseed
171
      DEcformals$initialpop <- seed
172
      rpseed <- FALSE
173
    } else {
174
      rpseed <- TRUE
175
    }
176
	if(hasArg(rpseed) & isTRUE(rpseed)) {
177
	    # initial seed population is generated with random_portfolios function
178
	    if(hasArg(eps)) eps=match.call(expand.dots=TRUE)$eps else eps = 0.01
179
    	rpconstraint<-constraint(assets=length(lower), min_sum=constraints$min_sum-eps, max_sum=constraints$max_sum+eps, 
180
             					min=lower, max=upper, weight_seq=generatesequence())
181
    	rp <- random_portfolios_v1(rpconstraints=rpconstraint,permutations=NP)
182
    	DEcformals$initialpop=rp
183
    }
184
    controlDE <- do.call(DEoptim::DEoptim.control,DEcformals)
185

186
    # minw = try(DEoptim( constrained_objective ,  lower = lower[1:N] , upper = upper[1:N] , control = controlDE, R=R, constraints=constraints, ...=...)) # add ,silent=TRUE here?
187
    minw = try(DEoptim::DEoptim( constrained_objective_v1 ,  lower = lower[1:N] , upper = upper[1:N] , control = controlDE, R=R, constraints=constraints, nargs = dotargs , ...=...)) # add ,silent=TRUE here?
188
 
189
    if(inherits(minw,"try-error")) { minw=NULL }
190
    if(is.null(minw)){
191
        message(paste("Optimizer was unable to find a solution for target"))
192
        return (paste("Optimizer was unable to find a solution for target"))
193
    }
194
    
195
    if(isTRUE(tmptrace)) trace <- tmptrace
196
    
197
    weights = as.vector( minw$optim$bestmem)
198
    weights <- normalize_weights(weights)
199
    names(weights) = colnames(R)
200

201
    out = list(weights=weights, objective_measures=constrained_objective_v1(w=weights,R=R,constraints,trace=TRUE)$objective_measures,out=minw$optim$bestval, call=call)
202
    if (isTRUE(trace)){
203
        out$DEoutput=minw
204
        out$DEoptim_objective_results<-try(get('.objectivestorage',envir=.storage),silent=TRUE)
205
        rm('.objectivestorage',envir=.storage)
206
    }
207
    
208
  } ## end case for DEoptim
209
  
210
  
211
  if(optimize_method=="random"){
212
      # call random_portfolios() with constraints and search_size to create matrix of portfolios
213
      if(missing(rp) | is.null(rp)){
214
          rp<-random_portfolios_v1(rpconstraints=constraints,permutations=search_size)
215
      }
216
      # store matrix in out if trace=TRUE
217
      if (isTRUE(trace)) out$random_portfolios<-rp
218
      # write foreach loop to call constrained_objective() with each portfolio
219
      if ("package:foreach" %in% search() & !hasArg(parallel)){
220
          ii <- 1
221
          rp_objective_results<-foreach::foreach(ii=1:nrow(rp), .errorhandling='pass') %dopar% constrained_objective_v1(w=rp[ii,],R,constraints,trace=trace,...=dotargs)
222
      } else {
223
          rp_objective_results<-apply(rp, 1, constrained_objective_v1, R=R, constraints=constraints, trace=trace, ...=dotargs)
224
      }
225
      # if trace=TRUE , store results of foreach in out$random_results
226
      if(isTRUE(trace)) out$random_portfolio_objective_results<-rp_objective_results
227
      # loop through results keeping track of the minimum value of rp_objective_results[[objective]]$out
228
      search<-vector(length=length(rp_objective_results))
229
      # first we construct the vector of results
230
      for (i in 1:length(search)) {
231
          if (isTRUE(trace)) {
232
              search[i]<-ifelse(try(rp_objective_results[[i]]$out),rp_objective_results[[i]]$out,1e6)
233
          } else {
234
              search[i]<-as.numeric(rp_objective_results[[i]])
235
          }
236
      }
237
      # now find the weights that correspond to the minimum score from the constrained objective
238
      # and normalize_weights so that we meet our min_sum/max_sum constraints
239
      if (isTRUE(trace)) {
240
          min_objective_weights<- try(normalize_weights(rp_objective_results[[which.min(search)]]$weights))
241
      } else {
242
          min_objective_weights<- try(normalize_weights(rp[which.min(search),]))
243
      }
244
      # re-call constrained_objective on the best portfolio, as above in DEoptim, with trace=TRUE to get results for out list
245
      out$weights<-min_objective_weights
246
      out$objective_measures<-try(constrained_objective_v1(w=min_objective_weights,R=R,constraints,trace=TRUE)$objective_measures)
247
      out$call<-call
248
      # construct out list to be as similar as possible to DEoptim list, within reason
249

250
  } ## end case for random
251
  
252
  
253
  if(optimize_method == "ROI_old"){
254
    # This will take a new constraint object that is of the same structure of a 
255
    # ROI constraint object, but with an additional solver arg.
256
    # then we can do something like this
257
    print("ROI_old is going to be depricated.")
258
    roi.result <- ROI::ROI_solve(x=constraints$constrainted_objective, constraints$solver)
259
    weights <- roi.result$solution
260
    names(weights) <- colnames(R)
261
    out$weights <- weights
262
    out$objective_measures <- roi.result$objval
263
    out$call <- call
264
  } ## end case for ROI_old
265
  
266
  
267
  
268
  if(optimize_method == "ROI"){
269
    # This takes in a regular constraint object and extracts the desired business objectives
270
    # and converts them to matrix form to be inputed into a closed form solver
271
    # Applying box constraints
272
    bnds <- list(lower = list(ind = seq.int(1L, N), val = as.numeric(constraints$min)),
273
                 upper = list(ind = seq.int(1L, N), val = as.numeric(constraints$max)))
274
    # retrive the objectives to minimize, these should either be "var" and/or "mean"
275
    # we can eight miniminze variance or maximize quiadratic utility (we will be minimizing the neg. quad. utility)
276
    moments <- list(mean=rep(0, N))
277
    alpha <- 0.05
278
    target <- NA
279
    for(objective in constraints$objectives){
280
      if(objective$enabled){
281
        if(!any(c(objective$name == "mean", objective$name == "var", objective$name == "CVaR")))
282
          stop("ROI only solves mean, var, or sample CVaR type business objectives, choose a different optimize_method.")
283
        # I'm not sure what changed, but moments$mean used to be a vector of the column means
284
        # now it is a scalar value of the mean of the entire R object
285
        if(objective$name == "mean"){
286
          moments[[objective$name]] <- try(as.vector(apply(R, 2, "mean", na.rm=TRUE)), silent=TRUE)
287
        } else {
288
          moments[[objective$name]] <- try(eval(as.symbol(objective$name))(R), silent=TRUE)
289
        }
290
        target <- ifelse(!is.null(objective$target),objective$target, target)
291
        alpha <- ifelse(!is.null(objective$alpha), objective$alpha, alpha)
292
        lambda <- ifelse(!is.null(objective$risk_aversion), objective$risk_aversion, 1)
293
      }
294
    }
295
    plugin <- ifelse(any(names(moments)=="var"), "quadprog", "glpk")  
296
    if(plugin == "quadprog") ROI_objective <- ROI::Q_objective(Q=2*lambda*moments$var, L=-moments$mean)
297
    if(plugin == "glpk") ROI_objective <- ROI::L_objective(L=-moments$mean)
298
    Amat <- rbind(rep(1, N), rep(1, N))
299
    dir.vec <- c(">=","<=")
300
    rhs.vec <- c(constraints$min_sum, constraints$max_sum)
301
    if(!is.na(target)) {
302
      Amat <- rbind(Amat, moments$mean)
303
      dir.vec <- c(dir.vec, "==")
304
      rhs.vec <- c(rhs.vec, target)
305
    }
306
    if(try(!is.null(constraints$groups), silent=TRUE)){
307
      if(sum(constraints$groups) != N)
308
        stop("Number of assets in each group needs to sum to number of total assets.")
309
      n.groups <- length(constraints$groups)
310
      if(!all(c(length(constraints$cLO),length(constraints$cLO)) == n.groups) )
311
        stop("Number of group constraints exceeds number of groups.")
312
      Amat.group <- matrix(0, nrow=n.groups, ncol=N)
313
      k <- 1
314
      l <- 0
315
      for(i in 1:n.groups){
316
        j <- constraints$groups[i] 
317
        Amat.group[i, k:(l+j)] <- 1
318
        k <- l + j + 1
319
        l <- k - 1
320
      }
321
      if(is.null(constraints$cLO)) cLO <- rep(-Inf, n.groups)
322
      if(is.null(constraints$cUP)) cUP <- rep(Inf, n.groups)
323
      Amat <- rbind(Amat, Amat.group, -Amat.group)
324
      dir.vec <- c(dir.vec, rep(">=", (n.groups + n.groups)))
325
      rhs.vec <- c(rhs.vec, constraints$cLO, -constraints$cUP)
326
    }
327
    if(any(names(moments)=="CVaR")) {
328
      Rmin <- ifelse(is.na(target), 0, target)
329
      ROI_objective <- ROI::L_objective(c(rep(0,N), rep(1/(alpha*T),T), 1))
330
      Amat <- cbind(rbind(1, 1, moments$mean, coredata(R)), rbind(0, 0, 0, cbind(diag(T), 1))) 
331
      dir.vec <- c(">=","<=",">=",rep(">=",T))
332
      rhs.vec <- c(constraints$min_sum, constraints$max_sum, Rmin ,rep(0, T))
333
      if(try(!is.null(constraints$groups), silent=TRUE)){
334
        zeros <- matrix(0, nrow=n.groups, ncol=(T+1))
335
        Amat <- rbind(Amat, cbind(Amat.group, zeros), cbind(-Amat.group, zeros))
336
        dir.vec <- c(dir.vec, rep(">=", (n.groups + n.groups)))
337
        rhs.vec <- c(rhs.vec, constraints$cLO, -constraints$cUP)
338
      }
339
    }
340
    opt.prob <- ROI::OP(objective=ROI_objective, 
341
                         constraints=ROI::L_constraint(L=Amat, dir=dir.vec, rhs=rhs.vec),
342
                         bounds=bnds)
343
    roi.result <- ROI::ROI_solve(x=opt.prob, solver=plugin)
344
    weights <- roi.result$solution[1:N]
345
    names(weights) <- colnames(R)
346
    out$weights <- weights
347
    out$out <- roi.result$objval
348
    out$call <- call
349
  } ## end case for ROI
350

351
  
352
  ## case if method=pso---particle swarm
353
  if(optimize_method=="pso"){
354
    stopifnot("package:pso" %in% search()  ||  requireNamespace("pso",quietly = TRUE) )
355
    if(hasArg(maxit)) maxit=match.call(expand.dots=TRUE)$maxit else maxit=N*50
356
    controlPSO <- list(trace=FALSE, fnscale=1, maxit=1000, maxf=Inf, abstol=-Inf, reltol=0)
357
    PSOcargs <- names(controlPSO)
358
    
359
    if( is.list(dotargs) ){
360
      pm <- pmatch(names(dotargs), PSOcargs, nomatch = 0L)
361
      names(dotargs[pm > 0L]) <- PSOcargs[pm]
362
      controlPSO$maxit <- maxit
363
      controlPSO[pm] <- dotargs[pm > 0L]
364
      if(!hasArg(reltol)) controlPSO$reltol <- .0001 # 1/100 of 1% change in objective is insignificant enough to restart a swarm
365
      #NOTE reltol has a different meaning for pso than it has for DEoptim.  for DEoptim, reltol is a stopping criteria, for pso, 
366
      #     it is a restart criteria.
367
        
368
      if(!hasArg(s)) {
369
        s <- N*10
370
        controlPSO$s<-s
371
        } #swarm size
372
      if(!hasArg(maxit.stagnate)) {
373
        #stopping criteria 
374
        maxit.stagnate <- controlPSO$s
375
        controlPSO$maxit.stagnate <- maxit.stagnate
376
        }     
377
      if(hasArg(trace) && try(trace==TRUE,silent=TRUE)) controlPSO$trace <- TRUE
378
      if(hasArg(trace) && isTRUE(trace)) {
379
          controlPSO$trace <- TRUE
380
          controlPSO$trace.stats=TRUE
381
      }
382
  }
383
    
384
    # get upper and lower weights parameters from constraints
385
    upper <- constraints$max
386
    lower <- constraints$min
387
    
388
    minw = try(pso::psoptim( par = rep(NA, N), fn = constrained_objective_v1 ,  R=R, constraints=constraints,
389
                        lower = lower[1:N] , upper = upper[1:N] , control = controlPSO)) # add ,silent=TRUE here?
390
    
391
    if(inherits(minw,"try-error")) { minw=NULL }
392
    if(is.null(minw)){
393
      message(paste("Optimizer was unable to find a solution for target"))
394
      return (paste("Optimizer was unable to find a solution for target"))
395
    }
396
    
397
    weights <- as.vector( minw$par)
398
    weights <- normalize_weights(weights)
399
    names(weights) <- colnames(R)
400
    
401
    out = list(weights=weights, 
402
               objective_measures=constrained_objective_v1(w=weights,R=R,constraints,trace=TRUE)$objective_measures,
403
               out=minw$value, 
404
               call=call)
405
    if (isTRUE(trace)){
406
      out$PSOoutput=minw
407
    }
408
    
409
  } ## end case for pso
410
  
411
  
412
  ## case if method=GenSA---Generalized Simulated Annealing
413
  if(optimize_method=="GenSA"){
414
    stopifnot("package:GenSA" %in% search()  ||  requireNamespace("GenSA",quietly = TRUE) )
415
    if(hasArg(maxit)) maxit=match.call(expand.dots=TRUE)$maxit else maxit=N*50
416
    controlGenSA <- list(maxit = 5000, threshold.stop = NULL, temperature = 5230, 
417
                          visiting.param = 2.62, acceptance.param = -5, max.time = NULL, 
418
                          nb.stop.improvement = 1e+06, smooth = TRUE, max.call = 1e+07, 
419
                          verbose = FALSE)
420
    GenSAcargs <- names(controlGenSA)
421
    
422
    if( is.list(dotargs) ){
423
      pm <- pmatch(names(dotargs), GenSAcargs, nomatch = 0L)
424
      names(dotargs[pm > 0L]) <- GenSAcargs[pm]
425
      controlGenSA$maxit <- maxit
426
      controlGenSA[pm] <- dotargs[pm > 0L]
427
      if(hasArg(trace) && try(trace==TRUE,silent=TRUE)) controlGenSA$verbose <- TRUE
428
    }
429
    
430
    upper <- constraints$max
431
    lower <- constraints$min
432
    
433
    if(!is.null(rp)) par = rp[,1] else par = rep(1/N, N)
434
    
435
    minw = try(GenSA::GenSA( par=par, lower = lower[1:N] , upper = upper[1:N], control = controlGenSA, 
436
                      fn = constrained_objective_v1 ,  R=R, constraints=constraints)) # add ,silent=TRUE here?
437
    
438
    if(inherits(minw,"try-error")) { minw=NULL }
439
    if(is.null(minw)){
440
      message(paste("Optimizer was unable to find a solution for target"))
441
      return (paste("Optimizer was unable to find a solution for target"))
442
    }
443
    
444
    weights <- as.vector(minw$par)
445
    weights <- normalize_weights(weights)
446
    names(weights) <- colnames(R)
447
    
448
    out = list(weights=weights, 
449
               objective_measures=constrained_objective_v1(w=weights,R=R,constraints,trace=TRUE)$objective_measures,
450
               out=minw$value, 
451
               call=call)
452
    if (isTRUE(trace)){
453
      out$GenSAoutput=minw
454
    }
455
    
456
  } ## end case for GenSA
457
  
458
  
459
    end_t<-Sys.time()
460
    # print(c("elapsed time:",round(end_t-start_t,2),":diff:",round(diff,2), ":stats: ", round(out$stats,4), ":targets:",out$targets))
461
    message(c("elapsed time:",end_t-start_t))
462
    out$constraints<-constraints
463
    out$data_summary<-list(first=first(R),last=last(R))
464
    out$elapsed_time<-end_t-start_t
465
    out$end_t<-as.character(Sys.time())
466
    class(out)<-c(paste("optimize.portfolio",optimize_method,sep='.'),"optimize.portfolio")
467
    return(out)
468
}
469

470
.onLoad <- function(lib, pkg) {
471
  if(!exists('.storage'))
472
    .storage <<- new.env()
473
}
474

475
#' Constrained optimization of portfolios
476
#' 
477
#' This function aims to provide a wrapper for constrained optimization of 
478
#' portfolios that specify constraints and objectives.
479
#' 
480
#' @details
481
#' This function currently supports DEoptim, random portfolios, pso, GenSA, ROI, osqp, Rglpk, mco, and CVXR solvers as back ends.
482
#' Additional back end contributions for Rmetrics, ghyp, etc. would be welcome.
483
#'
484
#' When using random portfolios, search_size is precisely that, how many 
485
#' portfolios to test.  You need to make sure to set your feasible weights 
486
#' in generatesequence to make sure you have search_size unique 
487
#' portfolios to test, typically by manipulating the 'by' parameter 
488
#' to select something smaller than .01 
489
#' (I often use .002, as .001 seems like overkill)
490
#' 
491
#' When using DE, search_size is decomposed into two other parameters 
492
#' which it interacts with, NP and itermax.
493
#' 
494
#' NP, the number of members in each population, is set to cap at 2000 in 
495
#' DEoptim, and by default is the number of parameters (assets/weights) * 10.
496
#' 
497
#' itermax, if not passed in dots, defaults to the number of parameters (assets/weights) * 50.
498
#' 
499
#' When using GenSA and want to set \code{verbose=TRUE}, instead use \code{trace}. 
500
#' 
501
#' If \code{optimize_method="ROI"} is specified, a default solver will be 
502
#' selected based on the optimization problem. The \code{glpk} solver is the
503
#' default solver for LP and MILP optimization problems. The \code{quadprog} 
504
#' solver is the default solver for QP optimization problems. For example,
505
#' \code{optimize_method = "quadprog"} can be specified and the optimization
506
#' problem will be solved via ROI using the quadprog solver.
507
#' 
508
#' The extension to ROI solves a limited type of convex optimization problems:
509
#' \itemize{
510
#' \item Maxmimize portfolio return subject leverage, box, group, position limit, target mean return, and/or factor exposure constraints on weights.
511
#' \item Minimize portfolio variance subject to leverage, box, group, turnover, and/or factor exposure constraints (otherwise known as global minimum variance portfolio).
512
#' \item Minimize portfolio variance subject to leverage, box, group, and/or factor exposure constraints and a desired portfolio return.
513
#' \item Maximize quadratic utility subject to leverage, box, group, target mean return, turnover, and/or factor exposure constraints and risk aversion parameter.
514
#' (The risk aversion parameter is passed into \code{optimize.portfolio} as an added argument to the \code{portfolio} object).
515
#' \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}. 
516
#' 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.
517
#' \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.
518
#' \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}. 
519
#' 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.
520
#' }
521
#' These problems also support a weight_concentration objective where concentration
522
#' of weights as measured by HHI is added as a penalty term to the quadratic objective.
523
#' 
524
#' Because these convex optimization problem are standardized, there is no need for a penalty term. 
525
#' The \code{multiplier} argument in \code{\link{add.objective}} passed into the complete constraint object are ignored by the ROI solver.
526
#'
527
#' If \code{optimize_method="CVXR"} is specified, a default solver will be selected based on the optimization problem.
528
#' The default solver for Quadratic Programming will be \code{OSQP}, 
529
#' and the default solver for Linear Problem and Second-Order Cone Programming will be \code{SCS}.
530
#' Specified CVXR solver can be given by using \code{optimize_method=c("CVXR", "CVXRsolver")}.
531
#' CVXR supports some commercial solvers, including CBC, CPLEX, GUROBI and MOSEK, and some open source solvers, including GLPK, GLPK_MI, OSQP, SCS and ECOS.
532
#' For example, \code{optimize_method = c("CVXR", "ECOS")} can be specified and the optimization problem will be solved via CVXR using the ECOS solver.
533
#' 
534
#' The extension to CVXR solves a limited type of convex optimization problems:
535
#' \itemize{
536
#' \item Maxmimize portfolio mean return subject leverage, box, group, and/or target mean return constraints
537
#' \item Minimize portfolio variance subject to leverage, box, group, and/or target mean return constraints (otherwise known as global minimum variance portfolio).
538
#' \item Maximize quadratic utility subject to leverage, box, group, and/or target mean return constraints and risk aversion parameter.
539
#' (The default risk aversion is 1, and specified risk aversion could be given by \code{risk_aversion = 1}.
540
#' The risk aversion parameter is passed into \code{optimize.portfolio} as an added argument to the \code{portfolio} object.)
541
#' \item Minimize portfolio ES/ETL/CVaR optimization subject to leverage, box, group, and/or target mean return constraints and tail probability parameter.
542
#' (The default tail probability is 0.05, and specified tail probability could be given by \code{arguments = list(p=0.95)}.
543
#' The tail probability parameter is passed into \code{optimize.portfolio} as an added argument to the \code{portfolio} object.)
544
#' \item Minimize portfolio CSM optimization subject to leverage, box, group, and/or target mean return constraints and tail probability parameter.
545
#' (The default tail probability is 0.05, and specified tail probability could be given by \code{arguments = list(p=0.95)}.
546
#' The tail probability parameter is passed into \code{optimize.portfolio} as an added argument to the \code{portfolio} object.)
547
#' \item Maximize portfolio mean return per unit standard deviation (i.e. the Sharpe Ratio) subject to leverage, box, group, and/or target mean return constraints.
548
#' It should be specified by \code{maxSR=TRUE} in \code{optimize.portfolio} with both mean and var/StdDev objectives.
549
#' Otherwise, the default action is to maximize quadratic utility.
550
#' \item Maximize portfolio mean return per unit ES (i.e. the ES ratio/STARR) subject to leverage, box, group, and/or target mean return constraints.
551
#' It could be specified by \code{maxSTARR=TRUE} or \code{ESratio=TRUE} in \code{optimize.portfolio} with both mean and ES objectives.
552
#' The default action is to maximize ES ratio. If \code{maxSTARR=FALSE} or \code{ESratio=FALSE} is given, the action will be minimizing ES.
553
#' \item Maximize portfolio mean return per unit CSM (i.e. the CSM ratio) subject to leverage, box, group, and/or target mean return constraints.
554
#' It could be specified by \code{CSMratio=TRUE} in \code{optimize.portfolio} with both mean and CSM objectives.
555
#' The default action is to maximize CSM ratio. If \code{CSMratio=FALSE} is given, the action will be minimizing CSM.
556
#' }
557
#' 
558
#' Because these convex optimization problem are standardized, there is no need for a penalty term. 
559
#' The \code{multiplier} argument in \code{\link{add.objective}} passed into the complete constraint object are ignored by the CVXR solver.
560
#'
561
#'
562
#' @note
563
#' An object of class \code{v1_constraint} can be passed in for the \code{constraints} argument.
564
#' The \code{v1_constraint} object was used in the previous 'v1' specification to specify the 
565
#' constraints and objectives for the optimization problem, see \code{\link{constraint}}. 
566
#' We will attempt to detect if the object passed into the constraints argument 
567
#' is a \code{v1_constraint} object and update to the 'v2' specification by adding the 
568
#' constraints and objectives to the \code{portfolio} object.
569
#'
570
#' @param R an xts, vector, matrix, data frame, timeSeries or zoo object of asset returns
571
#' @param portfolio an object of type "portfolio" specifying the constraints and objectives for the optimization
572
#' @param constraints default=NULL, a list of constraint objects. An object of class 'v1_constraint' can be passed in here.
573
#' @param objectives default=NULL, a list of objective objects.
574
#' @param optimize_method one of "DEoptim", "random", "ROI", "pso", "GenSA", "osqp", "Rglpk", "mco", "CVXR", or a vector to specify CVXR solver.
575
#' A solver of ROI or CVXR can also be specified and will be solved via ROI or CVXR. See details.
576
#' @param search_size integer, how many portfolios to test, default 20,000
577
#' @param trace TRUE/FALSE if TRUE will attempt to return additional information on the path or portfolios searched
578
#' @param \dots any other passthru parameters
579
#' @param rp matrix of random portfolio weights, default NULL, mostly for automated use by rebalancing optimization or repeated tests on same portfolios
580
#' @param momentFUN the name of a function to call to set portfolio moments, default \code{\link{set.portfolio.moments_v2}}
581
#' @param message TRUE/FALSE. The default is message=FALSE. Display messages if TRUE.
582
#' 
583
#' @return a list containing the following elements
584
#' \describe{
585
#'   \item{\code{weights}:}{ The optimal set weights.}
586
#'   \item{\code{objective_measures}:}{ A list containing the value of each objective corresponding to the optimal weights.}
587
#'   \item{\code{opt_values}:}{ A list containing the value of each objective corresponding to the optimal weights.}
588
#'   \item{\code{out}:}{ The output of the solver.}
589
#'   \item{\code{call}:}{ The function call.}
590
#'   \item{\code{portfolio}:}{ The portfolio object.}
591
#'   \item{\code{R}:}{ The asset returns.}
592
#'   \item{\code{data summary:}}{ The first row and last row of \code{R}.}
593
#'   \item{\code{elapsed_time:}}{ The amount of time that elapses while the optimization is run.}
594
#'   \item{\code{end_t:}}{ The date and time the optimization completed.}
595
#' }
596
#' When Trace=TRUE is specified, the following elements will be returned in 
597
#' addition to the elements above. The output depends on the optimization 
598
#' method and is specific to each solver. Refer to the documentation of the
599
#' desired solver for more information.
600
#' 
601
#' \code{optimize_method="random"}
602
#' \describe{
603
#'   \item{\code{random_portfolios}:}{ A matrix of the random portfolios.}
604
#'   \item{\code{random_portfolio_objective_results}:}{ A list of the following elements for each random portfolio.}
605
#'   \describe{
606
#'     \item{\code{out}:}{ The output value of the solver corresponding to the random portfolio weights.}
607
#'     \item{\code{weights}:}{ The weights of the random portfolio.}
608
#'     \item{\code{objective_measures}:}{ A list of each objective measure corresponding to the random portfolio weights.}
609
#'   }
610
#' }
611
#' 
612
#' \code{optimize_method="DEoptim"}
613
#' \describe{
614
#'   \item{\code{DEoutput:}}{ A list (of length 2) containing the following elements:}
615
#'   \itemize{
616
#'     \item \code{optim}
617
#'     \item \code{member}
618
#'   }
619
#'   \item{\code{DEoptim_objective_results}:}{ A list containing the following elements for each intermediate population.}
620
#'   \itemize{
621
#'     \item \code{out}: The output of the solver.
622
#'     \item \code{weights}: Population weights.
623
#'     \item \code{init_weights}: Initial population weights.
624
#'     \item \code{objective_measures}: A list of each objective measure corresponding to the weights
625
#'   }
626
#' }
627
#' 
628
#' \code{optimize_method="pso"}
629
#' \itemize{
630
#'   \item \code{PSOoutput}: A list containing the following elements:
631
#'   \itemize{
632
#'     \item par
633
#'     \item value 
634
#'     \item counts 
635
#'     \item convergence
636
#'     \item message
637
#'     \item stats
638
#'   }
639
#' }
640
#' 
641
#' \code{optimize_method="GenSA"}
642
#' \itemize{
643
#'   \item \code{GenSAoutput}: A list containing the following elements:
644
#'   \itemize{
645
#'     \item value
646
#'     \item par
647
#'     \item trace.mat
648
#'     \item counts
649
#'   }
650
#' }
651
#' 
652
#' @author Kris Boudt, Peter Carl, Brian G. Peterson, Ross Bennett, Xiaokang Feng, Xinran Zhao
653
#' @aliases optimize.portfolio_v2 optimize.portfolio_v1
654
#' @seealso \code{\link{portfolio.spec}}
655
#' @rdname optimize.portfolio
656
#' @name optimize.portfolio
657
#' @export optimize.portfolio
658
#' @export optimize.portfolio_v2
659
optimize.portfolio <- optimize.portfolio_v2 <- function(
660
  R,
661
  portfolio=NULL,
662
  constraints=NULL,
663
  objectives=NULL,
664
  optimize_method=c("DEoptim","random","ROI","pso","GenSA", "Rglpk", "osqp", "mco", "CVXR", "cvxr", ...),
665
  search_size=20000,
666
  trace=FALSE, ...,
667
  rp=NULL,
668
  momentFUN='set.portfolio.moments',
669
  message=FALSE
670
)
671
{
672
  # This is the case where the user has passed in a list of portfolio objects
673
  # for the portfolio argument.
674
  # Loop through the portfolio list and recursively call optimize.portfolio
675
  # Note that I return at the end of this block. I know it is not good practice
676
  # to return before the end of a function, but I am not sure of another way
677
  # to handle a list of portfolio objects with the recursive call to 
678
  # optimize.portfolio. 
679
  if(inherits(portfolio, "portfolio.list")){
680
    n.portf <- length(portfolio)
681
    opt.list <- vector("list", n.portf)
682
    for(i in 1:length(opt.list)){
683
      if(message) cat("Starting optimization of portfolio ", i, "\n")
684
      opt.list[[i]] <- optimize.portfolio(R=R, 
685
                                          portfolio=portfolio[[i]],
686
                                          constraints=constraints, 
687
                                          objectives=objectives, 
688
                                          optimize_method=optimize_method, 
689
                                          search_size=search_size, 
690
                                          trace=trace, 
691
                                          ...=..., 
692
                                          rp=rp, 
693
                                          momentFUN=momentFUN, 
694
                                          message=message)
695
    }
696
    out <- combine.optimizations(opt.list)
697
    ##### return here for portfolio.list because this is a recursive call
698
    ##### for optimize.portfolio
699
    return(out)
700
  }
701
  
702
  # Detect regime switching portfolio
703
  if(inherits(portfolio, "regime.portfolios")){
704
    regime.switching <- TRUE
705
    regime <- portfolio$regime
706
    if(index(last(R)) %in% index(regime)){
707
      regime.idx <- as.numeric(regime[index(last(R))])[1]
708
      portfolio <- portfolio$portfolio.list[[regime.idx]]
709
      #cat("regime: ", regime.idx, "\n")
710
    } else {
711
      warning("Dates in regime and R do not match, defaulting to first portfolio")
712
      regime.idx <- 1
713
      portfolio <- portfolio$portfolio.list[[regime.idx]]
714
    }
715
  } else {
716
    regime.switching <- FALSE
717
  }
718
  
719
  # This is the case where the user has passed in a mult.portfolio.spec
720
  # object for multiple layer portfolio optimization.
721
  if(inherits(portfolio, "mult.portfolio.spec")){
722
    # This function calls optimize.portfolio.rebalancing on each sub portfolio
723
    # according to the given optimization parameters and returns an xts object
724
    # representing the proxy returns of each sub portfolio.
725
    R <- proxy.mult.portfolio(R=R, mult.portfolio=portfolio)
726
    
727
    # The optimization is controlled by the constraints and objectives in the
728
    # top level portfolio so now set the 'portfolio' to the top level portfolio
729
    portfolio <- portfolio$top.portfolio
730
  }
731
  
732
  # Optimization Model Language
733
  if(length(optimize_method) == 2) optimize_method <- optimize_method[2] else optimize_method <- optimize_method[1]
734
  
735
  tmptrace <- NULL
736
  start_t <- Sys.time()
737
  
738
  #store the call for later
739
  call <- match.call()
740
  
741
  if (!is.null(portfolio) & !is.portfolio(portfolio)){
742
    stop("you must pass in an object of class 'portfolio' to control the optimization")
743
  }
744
  
745
  # Check for constraints and objectives passed in separately outside of the portfolio object
746
  if(!is.null(constraints)){
747
    if(inherits(constraints, "v1_constraint")){
748
      if(is.null(portfolio)){
749
        # If the user has not passed in a portfolio, we will create one for them
750
        tmp_portf <- portfolio.spec(assets=constraints$assets)
751
      }
752
      message("constraint object passed in is a 'v1_constraint' object, updating to v2 specification")
753
      portfolio <- update_constraint_v1tov2(portfolio=tmp_portf, v1_constraint=constraints)
754
      # print.default(portfolio)
755
    }
756
    if(!inherits(constraints, "v1_constraint")){
757
      # Insert the constraints into the portfolio object
758
      portfolio <- insert_constraints(portfolio=portfolio, constraints=constraints)
759
    }
760
  }
761
  if(!is.null(objectives)){
762
    # Insert the objectives into the portfolio object
763
    portfolio <- insert_objectives(portfolio=portfolio, objectives=objectives)
764
  }
765
  
766
  R <- checkData(R)
767
  N <- length(portfolio$assets)
768
  if (ncol(R) > N) {
769
    R <- R[,names(portfolio$assets)]
770
  }
771
  T <- nrow(R)
772
  
773
  # Initialize an empty list used as the return object
774
  out <- list()
775
  
776
  weights <- NULL 
777
  
778
  # Get the constraints from the portfolio object
779
  constraints <- get_constraints(portfolio)
780
  
781
  # set portfolio moments only once
782
  # For set.portfolio.moments, we are passing the returns,
783
  # portfolio object, and dotargs. dotargs is a list of arguments
784
  # that are passed in as dots in optimize.portfolio. This was
785
  # causing errors if clean="boudt" was specified in an objective
786
  # and an argument such as itermax was passed in as dots to 
787
  # optimize.portfolio. See r2931
788
  moment_name = momentFUN
789
  if(!is.function(momentFUN)){
790
    momentFUN <- match.fun(momentFUN)
791
  }
792
  
793
  # **
794
  # When an ES/ETL/CVaR problem is being solved by a linear solver, the higher
795
  # moments do not need to be calculated. The moments are very compute 
796
  # intensive and slow down the optimization problem.
797
  
798
  # match the args for momentFUN
799
  .formals <- formals(momentFUN)
800
  .formals <- modify.args(formals=.formals, arglist=list(...), dots=TRUE)
801
  # ** pass ROI=TRUE to set.portfolio.moments so the moments are not calculated
802
  if(optimize_method %in% c("ROI", "quadprog", "glpk", "symphony", "ipop")){
803
    obj_names <- unlist(lapply(portfolio$objectives, function(x) x$name))
804
    if(any(obj_names %in% c("CVaR", "ES", "ETL"))){
805
      .formals <- modify.args(formals=.formals, arglist=list(ROI=TRUE), dots=TRUE)
806
    }
807
  }
808
  if("R" %in% names(.formals)) .formals <- modify.args(formals=.formals, arglist=NULL, R=R, dots=FALSE)
809
  if("portfolio" %in% names(.formals)) .formals <- modify.args(formals=.formals, arglist=NULL, portfolio=portfolio, dots=FALSE)
810
  .formals$... <- NULL
811
  
812
  # call momentFUN
813
  mout <- try(do.call(momentFUN, .formals), silent=TRUE)
814
  
815
  if(inherits(mout, "try-error")) { 
816
    message(paste("portfolio moment function failed with message", mout))
817
  } else {
818
    #.args_env <- as.environment(mout)
819
    #.args_env <- new.env()
820
    # Assign each element of mout to the .args_env environment
821
    #for(name in names(mout)){
822
    #  .args_env[[name]] <- mout[[name]]
823
    #}
824
    dotargs <- mout
825
  }
826
  
827
  # Function to normalize weights to min_sum and max_sum
828
  # This function could be replaced by rp_transform
829
  normalize_weights <- function(weights){
830
    # normalize results if necessary
831
    if(!is.null(constraints$min_sum) | !is.null(constraints$max_sum)){
832
      # the user has passed in either min_sum or max_sum constraints for the portfolio, or both.
833
      # we'll normalize the weights passed in to whichever boundary condition has been violated
834
      # NOTE: this means that the weights produced by a numeric optimization algorithm like DEoptim
835
      # might violate your constraints, so you'd need to renormalize them after optimizing
836
      # we'll create functions for that so the user is less likely to mess it up.
837
      
838
      # NOTE: need to normalize in the optimization wrapper too before we return, since we've normalized in here
839
      # In Kris' original function, this was manifested as a full investment constraint
840
      if(!is.null(constraints$max_sum) & constraints$max_sum != Inf & constraints$max_sum != 0) {
841
        max_sum=constraints$max_sum
842
        if(sum(weights)>max_sum) { weights<-(max_sum/sum(weights))*weights } # normalize to max_sum
843
      }
844
      
845
      if(!is.null(constraints$min_sum) & constraints$min_sum != -Inf & constraints$min_sum != 0) {
846
        min_sum=constraints$min_sum
847
        if(sum(weights)<min_sum) { weights<-(min_sum/sum(weights))*weights } # normalize to min_sum
848
      }
849
      
850
    } # end min_sum and max_sum normalization
851
    return(weights)
852
  }
853
  
854
  # DEoptim optimization method
855
  if(optimize_method == "DEoptim"){
856
    stopifnot("package:DEoptim" %in% search()  ||  requireNamespace("DEoptim",quietly = TRUE))
857
    # DEoptim does 200 generations by default, so lets set the size of each generation to search_size/200)
858
    if(hasArg(itermax)) itermax=match.call(expand.dots=TRUE)$itermax else itermax=N*50
859
    NP <- round(search_size/itermax)
860
    if(NP < (N * 10)) NP <- N * 10
861
    if(NP > 2000) NP <- 2000
862
    if(!hasArg(itermax)) {
863
      itermax <- round(search_size / NP)
864
      if(itermax < 50) itermax <- 50 #set minimum number of generations
865
    }
866
    
867
    #check to see whether we need to disable foreach for parallel optimization, esp if called from inside foreach
868
    if(hasArg(parallel)) parallel <- match.call(expand.dots=TRUE)$parallel else parallel <- TRUE
869
    if(!isTRUE(parallel) && 'package:foreach' %in% search()){
870
      foreach::registerDoSEQ()
871
    }
872
    
873
    DEcformals <- formals(DEoptim::DEoptim.control)
874
    DEcargs <- names(DEcformals)
875
    if( is.list(dotargs) ){
876
      pm <- pmatch(names(dotargs), DEcargs, nomatch = 0L)
877
      names(dotargs[pm > 0L]) <- DEcargs[pm]
878
      DEcformals$NP <- NP
879
      DEcformals$itermax <- itermax
880
      DEcformals[pm] <- dotargs[pm > 0L]
881
      if(!hasArg(strategy)) {
882
        # use DE/current-to-p-best/1
883
        strategy=6
884
        DEcformals$strategy=strategy
885
        }
886
      if(!hasArg(reltol)) {
887
        # 1/1000 of 1% change in objective is significant
888
        reltol=0.000001
889
        DEcformals$reltol=reltol
890
        } 
891
      if(!hasArg(steptol)) {
892
        # number of assets times 1.5 tries to improve
893
        steptol=round(N*1.5)
894
        DEcformals$steptol=steptol
895
        } 
896
      if(!hasArg(c)) {
897
        # JADE mutation parameter, this could maybe use some adjustment
898
        tmp.c=0.4
899
        DEcformals$c=tmp.c
900
        }
901
      if(!hasArg(storepopfrom)) {
902
        storepopfrom=1
903
        DEcformals$storepopfrom=storepopfrom
904
      }
905
      if(isTRUE(parallel) && 'package:foreach' %in% search()){
906
        if(!hasArg(parallelType)) {
907
          #use all cores
908
          parallelType=2
909
          DEcformals$parallelType=parallelType
910
          }
911
        if(!hasArg(packages)) {
912
          #use all packages
913
          packages <- names(sessionInfo()$otherPkgs)
914
          DEcformals$packages <- packages
915
          }
916
      }
917
      #TODO FIXME also check for a passed in controlDE list, including checking its class, and match formals
918
    }
919
    if(hasArg(traceDE)) traceDE=match.call(expand.dots=TRUE)$traceDE else traceDE=TRUE
920
    DEcformals$trace <- traceDE
921
    if(isTRUE(trace)) { 
922
      #we can't pass trace=TRUE into constrained objective with DEoptim, because it expects a single numeric return
923
      tmptrace <- trace 
924
      assign('.objectivestorage', list(), envir=as.environment(.storage))
925
      trace=FALSE
926
    } 
927
    
928
    # get upper and lower weights parameters from constraints
929
    upper <- constraints$max
930
    lower <- constraints$min
931
    
932
    # issue message if min_sum and max_sum are restrictive
933
    if((constraints$max_sum - constraints$min_sum) < 0.02){
934
      message("Leverage constraint min_sum and max_sum are restrictive, 
935
              consider relaxing. e.g. 'full_investment' constraint should be min_sum=0.99 and max_sum=1.01")
936
    }
937
    
938
    #if(hasArg(rpseed)){ 
939
    #  seed <- match.call(expand.dots=TRUE)$rpseed
940
    #  DEcformals$initialpop <- seed
941
    #  rpseed <- FALSE
942
    #} else {
943
    #  rpseed <- TRUE
944
    #}
945
    #if(hasArg(rpseed) & isTRUE(rpseed)) {
946
    #  # initial seed population is generated with random_portfolios function
947
    #  # if(hasArg(eps)) eps=match.call(expand.dots=TRUE)$eps else eps = 0.01
948
    #  if(hasArg(rp_method)) rp_method=match.call(expand.dots=TRUE)$rp_method else rp_method="sample"
949
    #  if(hasArg(eliminate)) eliminate=match.call(expand.dots=TRUE)$eliminate else eliminate=TRUE
950
    #  if(hasArg(fev)) fev=match.call(expand.dots=TRUE)$fev else fev=0:5
951
    #  rp <- random_portfolios(portfolio=portfolio, permutations=NP, rp_method=rp_method, eliminate=eliminate, fev=fev)
952
    #  DEcformals$initialpop <- rp
953
    #}
954
    
955
    # Use rp as the initial population or generate from random portfolios
956
    if(!is.null(rp)){
957
      rp_len <- min(nrow(rp), NP)
958
      seed <- rp[1:rp_len,]
959
      DEcformals$initialpop <- seed
960
    } else{
961
      # Initial seed population is generated with random_portfolios function if rp is not passed in
962
      if(hasArg(rp_method)) rp_method=match.call(expand.dots=TRUE)$rp_method else rp_method="sample"
963
      # if(hasArg(eliminate)) eliminate=match.call(expand.dots=TRUE)$eliminate else eliminate=TRUE
964
      if(hasArg(fev)) fev=match.call(expand.dots=TRUE)$fev else fev=0:5
965
      rp <- random_portfolios(portfolio=portfolio, permutations=(NP+1), rp_method=rp_method, eliminate=FALSE, fev=fev)
966
      DEcformals$initialpop <- rp
967
    }
968
    
969
    controlDE <- do.call(DEoptim::DEoptim.control, DEcformals)
970
    # We are passing fn_map to the optional fnMap function to do the 
971
    # transformation so we need to force normalize=FALSE in call to constrained_objective
972
    minw = try(DEoptim::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)
973
    
974
    if(inherits(minw, "try-error")) { 
975
      message(minw)
976
      minw=NULL
977
    }
978
    if(is.null(minw)){
979
      message(paste("Optimizer was unable to find a solution for target"))
980
      return (paste("Optimizer was unable to find a solution for target"))
981
    }
982
    
983
    if(isTRUE(tmptrace)) trace <- tmptrace
984
    
985
    weights <- as.vector(minw$optim$bestmem)
986
    print(weights)
987
    # is it necessary to normalize the weights here?
988
    # weights <- normalize_weights(weights)
989
    names(weights) <- colnames(R)
990
    obj_vals <- constrained_objective(w=weights, R=R, portfolio, trace=TRUE, normalize=FALSE, env=dotargs)$objective_measures
991
    out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=minw$optim$bestval, call=call)
992
    if (isTRUE(trace)){
993
      out$DEoutput <- minw
994
      out$DEoptim_objective_results <- try(get('.objectivestorage',envir=.storage),silent=TRUE)
995
      rm('.objectivestorage',envir=.storage)
996
      #out$DEoptim_objective_results <- try(get('.objectivestorage',pos='.GlobalEnv'),silent=TRUE)
997
      #rm('.objectivestorage',pos='.GlobalEnv')
998
    }
999
    
1000
  } ## end case for DEoptim
1001
  
1002
  # case for random portfolios optimization method
1003
  if(optimize_method=="random"){
1004
    # issue message if min_sum and max_sum are too restrictive
1005
    if((constraints$max_sum - constraints$min_sum) < 0.02){
1006
      message("Leverage constraint min_sum and max_sum are restrictive, 
1007
              consider relaxing. e.g. 'full_investment' constraint should be min_sum=0.99 and max_sum=1.01")
1008
    }
1009
    
1010
    # call random_portfolios() with portfolio and search_size to create matrix of portfolios
1011
    if(missing(rp) | is.null(rp)){
1012
      if(hasArg(rp_method)) rp_method=match.call(expand.dots=TRUE)$rp_method else rp_method="sample"
1013
      if(hasArg(eliminate)) eliminate=match.call(expand.dots=TRUE)$eliminate else eliminate=TRUE
1014
      if(hasArg(fev)) fev=match.call(expand.dots=TRUE)$fev else fev=0:5
1015
      rp <- random_portfolios(portfolio=portfolio, permutations=search_size, rp_method=rp_method, eliminate=eliminate, fev=fev)
1016
    }
1017
    # store matrix in out if trace=TRUE
1018
    if (isTRUE(trace)) out$random_portfolios <- rp
1019
    # rp is already being generated with a call to fn_map so set normalize=FALSE in the call to constrained_objective
1020
    # write foreach loop to call constrained_objective() with each portfolio
1021
    if ("package:foreach" %in% search() & !hasArg(parallel)){
1022
      ii <- 1
1023
      rp_objective_results <- foreach::foreach(ii=1:nrow(rp), .errorhandling='pass') %dopar% constrained_objective(w=rp[ii,], R=R, portfolio=portfolio, trace=trace, env=dotargs, normalize=FALSE)
1024
    } else {
1025
      rp_objective_results <- apply(rp, 1, constrained_objective, R=R, portfolio=portfolio, trace=trace, normalize=FALSE, env=dotargs)
1026
    }
1027
    # if trace=TRUE , store results of foreach in out$random_results
1028
    if(isTRUE(trace)) out$random_portfolio_objective_results <- rp_objective_results
1029
    # loop through results keeping track of the minimum value of rp_objective_results[[objective]]$out
1030
    search <- vector(length=length(rp_objective_results))
1031
    # first we construct the vector of results
1032
    for (i in 1:length(search)) {
1033
      if (isTRUE(trace)) {
1034
        search[i] <- ifelse(try(rp_objective_results[[i]]$out), rp_objective_results[[i]]$out,1e6)
1035
      } else {
1036
        search[i] <- as.numeric(rp_objective_results[[i]])
1037
      }
1038
    }
1039
    # now find the weights that correspond to the minimum score from the constrained objective
1040
    # and normalize_weights so that we meet our min_sum/max_sum constraints
1041
    # Is it necessary to normalize the weights at all with random portfolios?
1042
    if (isTRUE(trace)) {
1043
      min_objective_weights <- try(normalize_weights(rp_objective_results[[which.min(search)]]$weights))
1044
    } else {
1045
      min_objective_weights <- try(normalize_weights(rp[which.min(search),]))
1046
    }
1047
    # re-call constrained_objective on the best portfolio, as above in DEoptim, with trace=TRUE to get results for out list
1048
    out$weights <- min_objective_weights
1049
    obj_vals <- try(constrained_objective(w=min_objective_weights, R=R, portfolio=portfolio, trace=TRUE, normalize=FALSE, env=dotargs)$objective_measures)
1050
    out$objective_measures <- obj_vals
1051
    out$opt_values <- obj_vals
1052
    out$call <- call
1053
    # construct out list to be as similar as possible to DEoptim list, within reason
1054
    
1055
  } ## end case for random
1056
  
1057
  roi_solvers <- c("ROI", "quadprog", "glpk", "symphony", "ipop")
1058
  if(optimize_method %in% roi_solvers){
1059
    # check for a control argument for list of solver control arguments
1060
    if(hasArg(control)) control=match.call(expand.dots=TRUE)$control else control=NULL
1061
    
1062
    # This takes in a regular portfolio object and extracts the constraints and
1063
    # objectives and converts them for input. to a closed form solver using
1064
    # ROI as an interface.
1065
    moments <- list(mean=rep(0, N))
1066
    alpha <- 0.05
1067
    if(!is.null(constraints$return_target)){
1068
      target <- constraints$return_target
1069
    } else {
1070
      target <- NA
1071
    }
1072
    lambda <- 1
1073
    
1074
    # list of valid objective names for ROI solvers
1075
    valid_objnames <- c("HHI", "mean", "var", "sd", "StdDev", "CVaR", "ES", "ETL")
1076
    #objnames <- unlist(lapply(portfolio$objectives, function(x) x$name))
1077
    for(objective in portfolio$objectives){
1078
      if(objective$enabled){
1079
        if(!(objective$name %in% valid_objnames)){
1080
          stop("ROI only solves mean, var/StdDev, HHI, or sample ETL/ES/CVaR type business objectives, choose a different optimize_method.")
1081
        }
1082
        
1083
        # Grab the arguments list per objective
1084
        # Currently we are only getting arguments for "p" and "clean", not sure if we need others for the ROI QP/LP solvers
1085
        # if(length(objective$arguments) >= 1) arguments <- objective$arguments else arguments <- list()
1086
        arguments <- objective$arguments
1087
        if(!is.null(arguments$clean)) clean <- arguments$clean else clean <- "none"
1088
        # Note: arguments$p grabs arguments$portfolio_method if no p is specified
1089
        # so we need to be explicit with arguments[["p"]]
1090
        if(!is.null(arguments[["p"]])) alpha <- arguments$p else alpha <- alpha
1091
        if(alpha > 0.5) alpha <- (1 - alpha)
1092
        
1093
        # Some of the sub-functions for optimizations use the returns object as
1094
        # part of the constraints matrix (e.g. etl_opt and etl_milp_opt) so we
1095
        # will store the cleaned returns in the moments object. This may not
1096
        # be the most efficient way to pass around a cleaned returns object, 
1097
        # but it will keep it separate from the R object passed in by the user
1098
        # and avoid "re-cleaning" already cleaned returns if specified in 
1099
        # multiple objectives.
1100
        if(clean != "none") moments$cleanR <- Return.clean(R=R, method=clean)
1101
        
1102
        # Use $mu and $sigma estimates from momentFUN if available, fall back to
1103
        # calculating sample mean and variance
1104
        if(objective$name == "mean"){
1105
          if(!is.null(mout$mu)){
1106
            moments[["mean"]] <- as.vector(mout$mu)
1107
          } else {
1108
            moments[["mean"]] <- try(as.vector(apply(Return.clean(R=R, method=clean), 2, "mean", na.rm=TRUE)), silent=TRUE)
1109
          }
1110
        } else if(objective$name %in% c("StdDev", "sd", "var")){
1111
          if(!is.null(mout$sigma)){
1112
            moments[["var"]] <- mout$sigma
1113
          } else {
1114
            moments[["var"]] <- try(var(x=Return.clean(R=R, method=clean), na.rm=TRUE), silent=TRUE)
1115
          }
1116
        } else if(objective$name %in% c("CVaR", "ES", "ETL")){
1117
          # do nothing (no point in computing ES here)
1118
          moments[[objective$name]] <- ""
1119
        } else {
1120
          # I'm not sure this is serving any purpose since we control the types
1121
          # of problems solved with ROI. The min ES problem only uses 
1122
          # moments$mu if a target return is specified.
1123
          moments[[objective$name]] <- try(eval(as.symbol(objective$name))(Return.clean(R=R, method=clean)), silent=TRUE)
1124
        }
1125
        target <- ifelse(!is.null(objective$target), objective$target, target)
1126
        # alpha <- ifelse(!is.null(objective$alpha), objective$alpha, alpha)
1127
        # only accept confidence level for ES/ETL/CVaR to come from the 
1128
        # arguments list to be consistent with how this is done in other solvers.
1129
        lambda <- ifelse(!is.null(objective$risk_aversion), objective$risk_aversion, lambda)
1130
        if(!is.null(objective$conc_aversion)) lambda_hhi <- objective$conc_aversion else lambda_hhi <- NULL
1131
        if(!is.null(objective$conc_groups)) conc_groups <- objective$conc_groups else conc_groups <- NULL
1132
      }
1133
    }
1134
    
1135
    if("var" %in% names(moments)){
1136
      # Set a default solver if optimize_method == "ROI", otherwise assume the
1137
      # optimize_method specified by the user is the solver for ROI
1138
      if(optimize_method == "ROI"){
1139
        solver <- "quadprog"
1140
      } else {
1141
        solver <- optimize_method
1142
      }
1143
      # Minimize variance if the only objective specified is variance
1144
      # Maximize Quadratic Utility if var and mean are specified as objectives
1145
      if(!is.null(constraints$turnover_target) | !is.null(constraints$ptc) | !is.null(constraints$leverage)){
1146
        if(!is.null(constraints$turnover_target) & !is.null(constraints$ptc)){
1147
          warning("Turnover and proportional transaction cost constraints detected, only running optimization for turnover constraint.")
1148
          constraints$ptc <- NULL
1149
        }
1150
        # turnover constraint
1151
        if(!is.null(constraints$turnover_target) & is.null(constraints$ptc)){
1152
          qp_result <- gmv_opt_toc(R=R, constraints=constraints, moments=moments, lambda=lambda, target=target, init_weights=portfolio$assets, solver=solver, control=control)
1153
          weights <- qp_result$weights
1154
          # obj_vals <- constrained_objective(w=weights, R=R, portfolio, trace=TRUE, normalize=FALSE)$objective_measures
1155
          obj_vals <- qp_result$obj_vals
1156
          out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=qp_result$out, call=call)
1157
        }
1158
        # proportional transaction costs constraint
1159
        if(!is.null(constraints$ptc) & is.null(constraints$turnover_target)){
1160
          qp_result <- gmv_opt_ptc(R=R, constraints=constraints, moments=moments, lambda=lambda, target=target, init_weights=portfolio$assets, solver=solver, control=control)
1161
          weights <- qp_result$weights
1162
          # obj_vals <- constrained_objective(w=weights, R=R, portfolio, trace=TRUE, normalize=FALSE)$objective_measures
1163
          obj_vals <- qp_result$obj_vals
1164
          out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=qp_result$out, call=call)
1165
        }
1166
        # leverage constraint
1167
        if(!is.null(constraints$leverage)){
1168
          qp_result <- gmv_opt_leverage(R=R, constraints=constraints, moments=moments, lambda=lambda, target=target, solver=solver, control=control)
1169
          weights <- qp_result$weights
1170
          # obj_vals <- constrained_objective(w=weights, R=R, portfolio, trace=TRUE, normalize=FALSE)$objective_measures
1171
          obj_vals <- qp_result$obj_vals
1172
          out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=qp_result$out, call=call)
1173
        }
1174
      } else {
1175
        # if(hasArg(ef)) ef=match.call(expand.dots=TRUE)$ef else ef=FALSE
1176
        if(hasArg(maxSR)) maxSR=match.call(expand.dots=TRUE)$maxSR else maxSR=FALSE
1177
        if(maxSR){
1178
          target <- max_sr_opt(R=R, constraints=constraints, moments=moments, lambda_hhi=lambda_hhi, conc_groups=conc_groups, solver=solver, control=control)
1179
          # need to set moments$mean=0 here because quadratic utility and target return is sensitive to returning no solution
1180
          tmp_moments_mean <- moments$mean
1181
          moments$mean <- rep(0, length(moments$mean))
1182
        }
1183
        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)
1184
        weights <- roi_result$weights
1185
        # obj_vals <- constrained_objective(w=weights, R=R, portfolio, trace=TRUE, normalize=FALSE)$objective_measures
1186
        obj_vals <- roi_result$obj_vals
1187
        if(maxSR){
1188
          # need to recalculate mean here if we are maximizing sharpe ratio
1189
          port.mean <- as.numeric(sum(weights * tmp_moments_mean))
1190
          names(port.mean) <- "mean"
1191
          obj_vals$mean <- port.mean
1192
        }
1193
        out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=roi_result$out, call=call)
1194
      }
1195
    }
1196
    if(length(names(moments)) == 1 & "mean" %in% names(moments)) {
1197
      # Set a default solver if optimize_method == "ROI", otherwise assume the
1198
      # optimize_method specified by the user is the solver for ROI
1199
      if(optimize_method == "ROI"){
1200
        solver <- "glpk"
1201
      } else {
1202
        solver <- optimize_method
1203
      }
1204
      
1205
      # Maximize return if the only objective specified is mean
1206
      if(!is.null(constraints$max_pos) | !is.null(constraints$leverage)) {
1207
        # This is an MILP problem if max_pos is specified as a constraint
1208
        roi_result <- maxret_milp_opt(R=R, constraints=constraints, moments=moments, target=target, solver=solver, control=control)
1209
        weights <- roi_result$weights
1210
        # obj_vals <- constrained_objective(w=weights, R=R, portfolio, trace=TRUE, normalize=FALSE)$objective_measures
1211
        obj_vals <- roi_result$obj_vals
1212
        out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=roi_result$out, call=call)
1213
      } else {
1214
        # Maximize return LP problem
1215
        roi_result <- maxret_opt(R=R, constraints=constraints, moments=moments, target=target, solver=solver, control=control)
1216
        weights <- roi_result$weights
1217
        # obj_vals <- constrained_objective(w=weights, R=R, portfolio, trace=TRUE, normalize=FALSE)$objective_measures
1218
        obj_vals <- roi_result$obj_vals
1219
        out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=roi_result$out, call=call)
1220
      }
1221
    }
1222
    if( any(c("CVaR", "ES", "ETL") %in% names(moments)) ) {
1223
      # Set a default solver if optimize_method == "ROI", otherwise assume the
1224
      # optimize_method specified by the user is the solver for ROI
1225
      if(optimize_method == "ROI"){
1226
        solver <- "glpk"
1227
      } else {
1228
        solver <- optimize_method
1229
      }
1230
      
1231
      if(hasArg(ef)) ef=match.call(expand.dots=TRUE)$ef else ef=FALSE
1232
      if(hasArg(maxSTARR)) maxSTARR=match.call(expand.dots=TRUE)$maxSTARR else maxSTARR=TRUE
1233
      if(ef) meanetl <- TRUE else meanetl <- FALSE
1234
      tmpnames <- c("CVaR", "ES", "ETL")
1235
      idx <- which(tmpnames %in% names(moments))
1236
      # Minimize sample ETL/ES/CVaR if CVaR, ETL, or ES is specified as an objective
1237
      if(length(moments) == 2 & all(moments$mean != 0) & ef==FALSE & maxSTARR){
1238
        # This is called by meanetl.efficient.frontier and we do not want that for efficient frontiers, need to have ef==FALSE
1239
        target <- mean_etl_opt(R=R, constraints=constraints, moments=moments, alpha=alpha, solver=solver, control=control)
1240
        meanetl <- TRUE
1241
      }
1242
      if(!is.null(constraints$max_pos)) {
1243
        # This is an MILP problem if max_pos is specified as a constraint
1244
        roi_result <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=target, alpha=alpha, solver=solver, control=control)
1245
        weights <- roi_result$weights
1246
        # obj_vals <- constrained_objective(w=weights, R=R, portfolio, trace=TRUE, normalize=FALSE)$objective_measures
1247
        # obj_vals <- roi_result$obj_vals
1248
        # calculate obj_vals based on solver output
1249
        obj_vals <- list()
1250
        if(meanetl) obj_vals$mean <- sum(weights * moments$mean)
1251
        obj_vals[[tmpnames[idx]]] <- roi_result$out
1252
        out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=roi_result$out, call=call)
1253
      } else {
1254
        # Minimize sample ETL/ES/CVaR LP Problem
1255
        roi_result <- etl_opt(R=R, constraints=constraints, moments=moments, target=target, alpha=alpha, solver=solver, control=control)
1256
        weights <- roi_result$weights
1257
        # obj_vals <- constrained_objective(w=weights, R=R, portfolio, trace=TRUE, normalize=FALSE)$objective_measures
1258
        # obj_vals <- roi_result$obj_vals
1259
        obj_vals <- list()
1260
        if(meanetl) obj_vals$mean <- sum(weights * moments$mean)
1261
        obj_vals[[tmpnames[idx]]] <- roi_result$out
1262
        out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=roi_result$out, call=call)
1263
      }
1264
    }
1265
    out$moment_values = list(momentFun = moment_name,mu = mout$mu, sigma = mout$sigma)
1266
    # Set here at the end so we get optimize.portfolio.ROI and not optimize.portfolio.{solver} classes
1267
    optimize_method <- "ROI"
1268
  } ## end case for ROI
1269
  
1270
  ## case if method=pso---particle swarm
1271
  if(optimize_method=="pso"){
1272
    stopifnot("package:pso" %in% search()  ||  requireNamespace("pso",quietly = TRUE) )
1273
    if(hasArg(maxit)) maxit=match.call(expand.dots=TRUE)$maxit else maxit=N*50
1274
    controlPSO <- list(trace=FALSE, fnscale=1, maxit=1000, maxf=Inf, abstol=-Inf, reltol=0)
1275
    PSOcargs <- names(controlPSO)
1276
    
1277
    if( is.list(dotargs) ){
1278
      pm <- pmatch(names(dotargs), PSOcargs, nomatch = 0L)
1279
      names(dotargs[pm > 0L]) <- PSOcargs[pm]
1280
      controlPSO$maxit <- maxit
1281
      controlPSO[pm] <- dotargs[pm > 0L]
1282
      if(!hasArg(reltol)) controlPSO$reltol <- .000001 # 1/1000 of 1% change in objective is significant
1283
      if(hasArg(trace) && try(trace==TRUE,silent=TRUE)) controlPSO$trace <- TRUE
1284
      if(hasArg(trace) && isTRUE(trace)) {
1285
        controlPSO$trace <- TRUE
1286
        controlPSO$trace.stats=TRUE
1287
      }
1288
    }
1289
    
1290
    # get upper and lower weights parameters from constraints
1291
    upper <- constraints$max
1292
    lower <- constraints$min
1293
    
1294
    minw <- try(pso::psoptim( par = rep(NA, N), fn = constrained_objective,  R=R, portfolio=portfolio, env=dotargs,
1295
                         lower = lower[1:N] , upper = upper[1:N] , control = controlPSO)) # add ,silent=TRUE here?
1296
    
1297
    if(inherits(minw,"try-error")) { minw=NULL }
1298
    if(is.null(minw)){
1299
      message(paste("Optimizer was unable to find a solution for target"))
1300
      return (paste("Optimizer was unable to find a solution for target"))
1301
    }
1302
    
1303
    weights <- as.vector( minw$par)
1304
    weights <- normalize_weights(weights)
1305
    names(weights) <- colnames(R)
1306
    obj_vals <- constrained_objective(w=weights, R=R, portfolio=portfolio, trace=TRUE, env=dotargs)$objective_measures
1307
    out <- list(weights=weights, 
1308
                objective_measures=obj_vals,
1309
                opt_values=obj_vals,
1310
                out=minw$value, 
1311
                call=call)
1312
    if (isTRUE(trace)){
1313
      out$PSOoutput=minw
1314
    }
1315
    
1316
  } ## end case for pso
1317
  
1318
  ## case if method=GenSA---Generalized Simulated Annealing
1319
  if(optimize_method=="GenSA"){
1320
    stopifnot("package:GenSA" %in% search()  ||  requireNamespace("GenSA",quietly = TRUE) )
1321
    if(hasArg(maxit)) maxit=match.call(expand.dots=TRUE)$maxit else maxit=N*50
1322
    controlGenSA <- list(maxit = 5000, threshold.stop = NULL, temperature = 5230, 
1323
                         visiting.param = 2.62, acceptance.param = -5, max.time = NULL, 
1324
                         nb.stop.improvement = 1e+06, smooth = TRUE, max.call = 1e+07, 
1325
                         verbose = FALSE)
1326
    GenSAcargs <- names(controlGenSA)
1327
    
1328
    if( is.list(dotargs) ){
1329
      pm <- pmatch(names(dotargs), GenSAcargs, nomatch = 0L)
1330
      names(dotargs[pm > 0L]) <- GenSAcargs[pm]
1331
      controlGenSA$maxit <- maxit
1332
      controlGenSA[pm] <- dotargs[pm > 0L]
1333
      if(hasArg(trace) && try(trace==TRUE,silent=TRUE)) controlGenSA$verbose <- TRUE
1334
    }
1335
    
1336
    upper <- constraints$max
1337
    lower <- constraints$min
1338

1339
    if(!is.null(rp)) par = rp[,1] else par = rep(1/N, N)
1340
      
1341
    minw = try(GenSA::GenSA( par=par, lower = lower[1:N] , upper = upper[1:N], control = controlGenSA, 
1342
                      fn = constrained_objective ,  R=R, portfolio=portfolio, env=dotargs)) # add ,silent=TRUE here?
1343
    
1344
    if(inherits(minw,"try-error")) { minw=NULL }
1345
    if(is.null(minw)){
1346
      message(paste("Optimizer was unable to find a solution for target"))
1347
      return (paste("Optimizer was unable to find a solution for target"))
1348
    }
1349
    
1350
    weights <- as.vector(minw$par)
1351
    weights <- normalize_weights(weights)
1352
    names(weights) <- colnames(R)
1353
    obj_vals <- constrained_objective(w=weights, R=R, portfolio=portfolio, trace=TRUE, env=dotargs)$objective_measures
1354
    out = list(weights=weights, 
1355
               objective_measures=obj_vals,
1356
               opt_values=obj_vals,
1357
               out=minw$value, 
1358
               call=call)
1359
    if (isTRUE(trace)){
1360
      out$GenSAoutput=minw
1361
    }
1362
    
1363
  } ## end case for GenSA
1364
  
1365
  ## case if method = Rglpk -- R GUI Linear Programming Kit
1366
  if (optimize_method == "Rglpk") {
1367
    # search for required package
1368
    stopifnot("package:Rglpk" %in% search() ||
1369
      requireNamespace("Rglpk", quietly = TRUE))
1370
  
1371
    # Rglpk solver can only solve linear programming problems
1372
    valid_objnames <- c("mean", "CVaR", "ES", "ETL")
1373
  
1374
    # default setting
1375
    target <- -Inf
1376
    alpha <- 0.05
1377
    reward <- FALSE
1378
    risk <- FALSE
1379
  
1380
    # search invalid objectives
1381
    for (objective in portfolio$objectives) {
1382
      if (objective$enabled) {
1383
        if (!(objective$name %in% valid_objnames)) {
1384
          stop("Rglpk only solves mean and ETL/ES/CVaR type business objectives, 
1385
                 choose a different optimize_method.")
1386
        }
1387
  
1388
        # alpha
1389
        alpha <- ifelse(!is.null(objective$arguments[["p"]]),
1390
          objective$arguments[["p"]], alpha
1391
        )
1392
  
1393
        # return target
1394
        target <- ifelse(!is.null(objective$target), objective$target, target)
1395
  
1396
        # optimization objective function
1397
        reward <- ifelse(objective$name == "mean", TRUE, reward)
1398
        risk <- ifelse(objective$name %in% valid_objnames[2:4], TRUE, risk)
1399
        arguments <- objective$arguments
1400
      }
1401
    }
1402
  
1403
    # trim alpha and target
1404
    alpha <- ifelse(alpha > 0.5, 1 - alpha, alpha)
1405
    target <- max(target, constraints$return_target, na.rm = TRUE)
1406
  
1407
    # get leverage paramters from constraints
1408
    max_sum <- constraints$max_sum
1409
    min_sum <- constraints$min_sum
1410
  
1411
    # transfer inf to big number
1412
    max_sum <- ifelse(is.infinite(max_sum), 9999.0, max_sum)
1413
    min_sum <- ifelse(is.infinite(min_sum), -9999.0, min_sum)
1414
  
1415
    # get upper and lower limitation
1416
    upper <- constraints$max
1417
    lower <- constraints$min
1418
  
1419
    upper[which(is.infinite(upper))] <- 9999.0
1420
    lower[which(is.infinite(lower))] <- 9999.0
1421
  
1422
    # issue message if min_sum and max_sum are restrictive
1423
    if ((max_sum - min_sum) < 0.02) {
1424
      message("Leverage constraint min_sum and max_sum are restrictive, 
1425
                consider relaxing. e.g. 'full_investment' constraint 
1426
                should be min_sum=0.99 and max_sum=1.01")
1427
    }
1428
  
1429
    # Here, I created two version optimization. One for no position
1430
    # limitation; the other one for situation with position limitation.
1431
    # This will keep code clean and accelerate simple case process speed.
1432
  
1433
    if (is.null(constraints$max_pos) & is.null(constraints$group_pos)) {
1434
      # deploy optimization for different types situations
1435
      # max return case
1436
      if (reward & !risk) {
1437
        # utility function coef
1438
        Rglpk.obj <- colMeans(R)
1439
  
1440
        # leverage constraint
1441
        Rglpk.mat <- rbind(rep(1, N), rep(1, N))
1442
        Rglpk.dir <- c(">=", "<=")
1443
        Rglpk.rhs <- c(min_sum, max_sum)
1444
  
1445
        # box constraint
1446
        Rglpk.bound <- list(
1447
          lower = list(
1448
            ind = 1:N,
1449
            val = lower
1450
          ),
1451
          upper = list(
1452
            ind = 1:N,
1453
            val = upper
1454
          )
1455
        )
1456
  
1457
        # group constraint
1458
        if (!is.null(constraints$groups)) {
1459
          # check group constraint validility
1460
          if (!all(c(
1461
            length(constraints$cLO),
1462
            length(constraints$cUP)
1463
          )
1464
          == length(constraints$groups))) {
1465
            stop("Please assign group constraint")
1466
          }
1467
  
1468
          # update constraints matrix, direction and right-hand vector
1469
          for (i in 1:length(constraints$groups)) {
1470
            # temp index variable
1471
            temp <- rep(0, N)
1472
            temp[constraints$groups[[i]]] <- 1
1473
  
1474
            Rglpk.mat <- rbind(Rglpk.mat, temp, temp)
1475
            Rglpk.dir <- c(Rglpk.dir, ">=", "<=")
1476
            Rglpk.rhs <- c(
1477
              Rglpk.rhs,
1478
              constraints$cLO[i],
1479
              constraints$cUP[i]
1480
            )
1481
          }
1482
          rm(temp)
1483
        }
1484
  
1485
        # return target constraint
1486
        if (!is.infinite(target)) {
1487
          Rglpk.mat <- rbind(Rglpk.mat, colMeans(R))
1488
          Rglpk.dir <- c(Rglpk.dir, ">=")
1489
          Rglpk.rhs <- c(Rglpk.rhs, target)
1490
        }
1491
  
1492
        # result from solver
1493
        Rglpk.result <- try(Rglpk::Rglpk_solve_LP(
1494
          obj = Rglpk.obj,
1495
          mat = Rglpk.mat,
1496
          dir = Rglpk.dir,
1497
          rhs = Rglpk.rhs,
1498
          bound = Rglpk.bound,
1499
          types = rep("C", N),
1500
          max = TRUE
1501
        ))
1502
  
1503
        # null result
1504
        if (inherits(Rglpk.result, "try-error")) Rglpk.result <- NULL
1505
        if (is.null(Rglpk.result)) {
1506
          message(paste("Optimizer was unable to find a solution for target"))
1507
          return(paste("Optimizer was unable to find a solution for target"))
1508
        }
1509
  
1510
        # prepare output
1511
        weights <- as.vector(Rglpk.result$solution[1:N])
1512
        weights <- normalize_weights(weights)
1513
        names(weights) <- colnames(R)
1514
        obj_vals <- constrained_objective(
1515
          w = weights, R = R, portfolio = portfolio,
1516
          trace = TRUE, env = dotargs
1517
        )$objective_measures
1518
        out <- list(
1519
          weights = weights,
1520
          objective_measures = obj_vals,
1521
          opt_values = obj_vals,
1522
          out = Rglpk.result$optimum,
1523
          call = call
1524
        )
1525
        if (isTRUE(trace)) {
1526
          out$Rglpkoutput <- Rglpk.result
1527
        }
1528
      }
1529
  
1530
      # min CVaR case
1531
      else if (risk & !reward) {
1532
        # utility function coef: weight + loss + VaR
1533
        Rglpk.obj <- c(
1534
          rep(0, N),
1535
          rep(-1 / alpha / T, T),
1536
          -1
1537
        )
1538
  
1539
        # leverage constraint
1540
        Rglpk.mat <- rbind(
1541
          c(rep(1, N), rep(0, T + 1)),
1542
          c(rep(1, N), rep(0, T + 1))
1543
        )
1544
        Rglpk.dir <- c(">=", "<=")
1545
        Rglpk.rhs <- c(min_sum, max_sum)
1546
  
1547
        # box constraint
1548
        Rglpk.bound <- list(
1549
          lower = list(
1550
            ind = c(1:N, N + T + 1),
1551
            val = c(lower, -Inf)
1552
          ),
1553
          upper = list(
1554
            ind = c(1:N, N + T + 1),
1555
            val = c(upper, Inf)
1556
          )
1557
        )
1558
  
1559
        # group constraint
1560
        if (!is.null(constraints$groups)) {
1561
          # check group constraint validility
1562
          if (!all(c(
1563
            length(constraints$cLO),
1564
            length(constraints$cUP)
1565
          )
1566
          == length(constraints$groups))) {
1567
            stop("Please assign group constraint")
1568
          }
1569
  
1570
          # update constraints matrix, direction and right-hand vector
1571
          for (i in 1:length(constraints$groups)) {
1572
            # temp index variable
1573
            temp <- rep(0, N + T + 1)
1574
            temp[constraints$groups[[i]]] <- 1
1575
  
1576
            Rglpk.mat <- rbind(Rglpk.mat, temp, temp)
1577
            Rglpk.dir <- c(Rglpk.dir, ">=", "<=")
1578
            Rglpk.rhs <- c(
1579
              Rglpk.rhs,
1580
              constraints$cLO[i],
1581
              constraints$cUP[i]
1582
            )
1583
          }
1584
          rm(temp)
1585
        }
1586
  
1587
        # return target constraint
1588
        if (!is.infinite(target)) {
1589
          Rglpk.mat <- rbind(Rglpk.mat, colMeans(R))
1590
          Rglpk.dir <- c(Rglpk.dir, ">=")
1591
          Rglpk.rhs <- c(Rglpk.rhs, target)
1592
        }
1593
  
1594
        # scenario constraint
1595
        Rglpk.mat <- rbind(
1596
          Rglpk.mat,
1597
          cbind(
1598
            matrix(R, T),
1599
            diag(1, T),
1600
            rep(1, T)
1601
          )
1602
        )
1603
        Rglpk.dir <- c(Rglpk.dir, rep(">=", T))
1604
        Rglpk.rhs <- c(Rglpk.rhs, rep(0, T))
1605
  
1606
        # result from solver
1607
        Rglpk.result <- try(Rglpk::Rglpk_solve_LP(
1608
          obj = Rglpk.obj,
1609
          mat = Rglpk.mat,
1610
          dir = Rglpk.dir,
1611
          rhs = Rglpk.rhs,
1612
          bound = Rglpk.bound,
1613
          types = rep("C", N + T + 1),
1614
          max = TRUE
1615
        ))
1616
  
1617
        # null result
1618
        if (inherits(Rglpk.result, "try-error")) Rglpk.result <- NULL
1619
        if (is.null(Rglpk.result)) {
1620
          message(paste("Optimizer was unable to find a solution for target"))
1621
          return(paste("Optimizer was unable to find a solution for target"))
1622
        }
1623
  
1624
        # prepare output
1625
        weights <- as.vector(Rglpk.result$solution[1:N])
1626
        weights <- normalize_weights(weights)
1627
        names(weights) <- colnames(R)
1628
        obj_vals <- constrained_objective(
1629
          w = weights, R = R, portfolio = portfolio,
1630
          trace = TRUE, env = dotargs
1631
        )$objective_measures
1632
        out <- list(
1633
          weights = weights,
1634
          objective_measures = obj_vals,
1635
          opt_values = obj_vals,
1636
          out = Rglpk.result$optimum,
1637
          call = call
1638
        )
1639
        if (isTRUE(trace)) {
1640
          out$Rglpkoutput <- Rglpk.result
1641
        }
1642
      }
1643
  
1644
      # max ratio
1645
      else if (risk & reward) {
1646
        # utility function coef: weight + loss + VaR + shrinkage
1647
        Rglpk.obj <- c(
1648
          rep(0, N), # weight
1649
          rep(-1 / alpha / T, T), # loss
1650
          -1, # VaR
1651
          0 # shrinkage
1652
        )
1653
  
1654
        # leverage constraint
1655
        Rglpk.mat <- rbind(
1656
          c(rep(1, N), rep(0, T + 1), -min_sum),
1657
          c(rep(1, N), rep(0, T + 1), -max_sum)
1658
        )
1659
        Rglpk.dir <- c(">=", "<=")
1660
        Rglpk.rhs <- c(0, 0)
1661
  
1662
        # box constraint
1663
        Rglpk.mat <- rbind(
1664
          Rglpk.mat,
1665
          cbind(diag(1, N), matrix(0, N, T + 1), -lower),
1666
          cbind(diag(1, N), matrix(0, N, T + 1), -upper)
1667
        )
1668
        Rglpk.dir <- c(Rglpk.dir, rep(">=", N), rep("<=", N))
1669
        Rglpk.rhs <- c(Rglpk.rhs, rep(0, 2 * N))
1670
  
1671
        Rglpk.bound <- list(
1672
          lower = list(
1673
            ind = c(1:N, N + T + 1),
1674
            val = c(rep(-Inf, N), -Inf)
1675
          ),
1676
          upper = list(
1677
            ind = c(1:N, N + T + 1),
1678
            val = c(rep(Inf, N), Inf)
1679
          )
1680
        )
1681
  
1682
        # group constraint
1683
        if (!is.null(constraints$groups)) {
1684
          # check group constraint validility
1685
          if (!all(c(
1686
            length(constraints$cLO),
1687
            length(constraints$cUP)
1688
          )
1689
          == length(constraints$groups))) {
1690
            stop("Please assign group constraint")
1691
          }
1692
  
1693
          # update constraints matrix, direction and right-hand vector
1694
          for (i in 1:length(constraints$groups)) {
1695
            # temp index variable
1696
            # low level
1697
            temp <- c(rep(0, N + T + 1), -constraints$cLO[i])
1698
            temp[constraints$groups[[i]]] <- 1
1699
            Rglpk.mat <- rbind(Rglpk.mat, temp)
1700
            Rglpk.dir <- c(Rglpk.dir, ">=")
1701
            Rglpk.rhs <- c(Rglpk.rhs, 0)
1702
  
1703
            # up level
1704
            temp <- c(rep(0, N + T + 1), -constraints$cUP[i])
1705
            temp[constraints$groups[[i]]] <- 1
1706
            Rglpk.mat <- rbind(Rglpk.mat, temp)
1707
            Rglpk.dir <- c(Rglpk.dir, "<=")
1708
            Rglpk.rhs <- c(Rglpk.rhs, 0)
1709
          }
1710
          rm(temp)
1711
        }
1712
  
1713
        # return target constraint
1714
        if (!is.infinite(target)) {
1715
          Rglpk.mat <- rbind(
1716
            Rglpk.mat,
1717
            c(rep(0, N + T + 1)),
1718
            target
1719
          )
1720
  
1721
          Rglpk.dir <- c(Rglpk.dir, "<=")
1722
          Rglpk.rhs <- c(Rglpk.rhs, 1)
1723
        }
1724
  
1725
        # shrinkage constraint
1726
        Rglpk.mat <- rbind(
1727
          Rglpk.mat,
1728
          c(colMeans(R), rep(0, T + 2))
1729
        )
1730
  
1731
        Rglpk.dir <- c(Rglpk.dir, "==")
1732
        Rglpk.rhs <- c(Rglpk.rhs, 1)
1733
  
1734
        # scenario constraint
1735
        Rglpk.mat <- rbind(
1736
          Rglpk.mat,
1737
          cbind(
1738
            matrix(R, T),
1739
            diag(1, T),
1740
            rep(1, T),
1741
            rep(0, T)
1742
          )
1743
        )
1744
        Rglpk.dir <- c(Rglpk.dir, rep(">=", T))
1745
        Rglpk.rhs <- c(Rglpk.rhs, rep(0, T))
1746
  
1747
  
1748
  
1749
        # result from solver
1750
        Rglpk.result <- try(Rglpk::Rglpk_solve_LP(
1751
          obj = Rglpk.obj,
1752
          mat = Rglpk.mat,
1753
          dir = Rglpk.dir,
1754
          rhs = Rglpk.rhs,
1755
          bound = Rglpk.bound,
1756
          types = rep("C", N + T + 2),
1757
          max = TRUE
1758
        ))
1759
  
1760
        # null result
1761
        if (inherits(Rglpk.result, "try-error")) Rglpk.result <- NULL
1762
        if (is.null(Rglpk.result)) {
1763
          message(paste("Optimizer was unable to find a solution for target"))
1764
          return(paste("Optimizer was unable to find a solution for target"))
1765
        }
1766
  
1767
        # prepare output
1768
        weights <- as.vector(Rglpk.result$solution[1:N] / Rglpk.result$solution[N + T + 2])
1769
        weights <- normalize_weights(weights)
1770
        names(weights) <- colnames(R)
1771
        obj_vals <- constrained_objective(
1772
          w = weights, R = R, portfolio = portfolio,
1773
          trace = TRUE, env = dotargs
1774
        )$objective_measures
1775
        out <- list(
1776
          weights = weights,
1777
          objective_measures = obj_vals,
1778
          opt_values = obj_vals,
1779
          out = Rglpk.result$optimum,
1780
          call = call
1781
        )
1782
        if (isTRUE(trace)) {
1783
          out$Rglpkoutput <- Rglpk.result
1784
        }
1785
      }
1786
    }
1787
  
1788
    # The simple case without position case is done.
1789
    # Following case can handle position limitation, but in a really slow
1790
    # pace, especially choosing n from 2n.
1791
  
1792
    else {
1793
      # deploy optimization for different types situations
1794
      # max return case
1795
      if (reward & !risk) {
1796
        # utility function coef: weight + position index
1797
        Rglpk.obj <- c(colMeans(R), rep(0, N))
1798
  
1799
        # leverage constraint
1800
        Rglpk.mat <- rbind(
1801
          c(rep(1, N), rep(0, N)),
1802
          c(rep(1, N), rep(0, N))
1803
        )
1804
        Rglpk.dir <- c(">=", "<=")
1805
        Rglpk.rhs <- c(min_sum, max_sum)
1806
  
1807
        # box constraint
1808
        Rglpk.bound <- list(
1809
          lower = list(
1810
            ind = 1:N,
1811
            val = lower
1812
          ),
1813
          upper = list(
1814
            ind = 1:N,
1815
            val = upper
1816
          )
1817
        )
1818
  
1819
        # group constraint
1820
        if (!is.null(constraints$groups)) {
1821
          # check group constraint validility
1822
          if (!all(c(
1823
            length(constraints$cLO),
1824
            length(constraints$cUP)
1825
          )
1826
          == length(constraints$groups))) {
1827
            stop("Please assign group constraint")
1828
          }
1829
  
1830
          # group weight constraint
1831
          for (i in 1:length(constraints$groups)) {
1832
            # temp index variable
1833
            temp <- rep(0, 2 * N)
1834
            temp[constraints$groups[[i]]] <- 1
1835
  
1836
            # weight constraints
1837
            Rglpk.mat <- rbind(Rglpk.mat, temp, temp)
1838
            Rglpk.dir <- c(Rglpk.dir, ">=", "<=")
1839
            Rglpk.rhs <- c(
1840
              Rglpk.rhs,
1841
              constraints$cLO[i],
1842
              constraints$cUP[i]
1843
            )
1844
          }
1845
  
1846
          # group postion limitation
1847
          if (!is.null(constraints$group_pos)) {
1848
            # check group position limitation length
1849
            if (length(constraints$group_pos)
1850
            != length(constraints$groups)) {
1851
              stop("Please assign group constraint")
1852
            }
1853
            for (i in 1:length(constraints$groups)) {
1854
              # temp index variable
1855
              temp <- rep(0, 2 * N)
1856
              temp[constraints$groups[[i]] + N] <- 1
1857
  
1858
              # position limitation constraints
1859
              Rglpk.mat <- rbind(Rglpk.mat, temp)
1860
              Rglpk.dir <- c(Rglpk.dir, "==")
1861
              Rglpk.rhs <- c(
1862
                Rglpk.rhs,
1863
                constraints$group_pos[i]
1864
              )
1865
            }
1866
          }
1867
  
1868
          # remove temp variable
1869
          rm(temp)
1870
        }
1871
  
1872
        # return target constraint
1873
        if (!is.infinite(target)) {
1874
          Rglpk.mat <- rbind(Rglpk.mat, c(colMeans(R), rep(0, N)))
1875
          Rglpk.dir <- c(Rglpk.dir, ">=")
1876
          Rglpk.rhs <- c(Rglpk.rhs, target)
1877
        }
1878
  
1879
        # position limitation constraint
1880
        if (!is.null(constraints$max_pos)) {
1881
          Rglpk.mat <- rbind(Rglpk.mat, c(rep(0, N), rep(1, N)))
1882
          Rglpk.dir <- c(Rglpk.dir, "<=")
1883
          Rglpk.rhs <- c(Rglpk.rhs, constraints$max_pos)
1884
  
1885
          Rglpk.mat <- rbind(
1886
            Rglpk.mat,
1887
            cbind(diag(1, N), diag(-upper, N)),
1888
            cbind(diag(1, N), diag(-lower, N))
1889
          )
1890
          Rglpk.dir <- c(Rglpk.dir, rep("<=", N), rep(">=", N))
1891
          Rglpk.rhs <- c(Rglpk.rhs, rep(0, 2 * N))
1892
        }
1893
  
1894
        # result from solver
1895
        Rglpk.result <- try(Rglpk::Rglpk_solve_LP(
1896
          obj = Rglpk.obj,
1897
          mat = Rglpk.mat,
1898
          dir = Rglpk.dir,
1899
          rhs = Rglpk.rhs,
1900
          bound = Rglpk.bound,
1901
          types = c(rep("C", N), rep("B", N)),
1902
          max = TRUE
1903
        ))
1904
  
1905
        # null result
1906
        if (inherits(Rglpk.result, "try-error")) Rglpk.result <- NULL
1907
        if (is.null(Rglpk.result)) {
1908
          message(paste("Optimizer was unable to find a solution for target"))
1909
          return(paste("Optimizer was unable to find a solution for target"))
1910
        }
1911
  
1912
        # prepare output
1913
        weights <- as.vector(Rglpk.result$solution[1:N])
1914
        weights <- normalize_weights(weights)
1915
        names(weights) <- colnames(R)
1916
        obj_vals <- constrained_objective(
1917
          w = weights, R = R, portfolio = portfolio,
1918
          trace = TRUE, env = dotargs
1919
        )$objective_measures
1920
        out <- list(
1921
          weights = weights,
1922
          objective_measures = obj_vals,
1923
          opt_values = obj_vals,
1924
          out = Rglpk.result$optimum,
1925
          call = call
1926
        )
1927
        if (isTRUE(trace)) {
1928
          out$Rglpkoutput <- Rglpk.result
1929
        }
1930
      }
1931
  
1932
      # min VaR case
1933
      else if (risk & !reward) {
1934
        # utility function coef: weight + loss + VaR + position index
1935
        Rglpk.obj <- c(
1936
          rep(0, N), # weight
1937
          rep(-1 / alpha / T, T), # loss
1938
          -1, # VaR
1939
          rep(0, N) # position index
1940
        )
1941
  
1942
        # leverage constraint
1943
        Rglpk.mat <- rbind(
1944
          c(rep(1, N), rep(0, T + 1 + N)),
1945
          c(rep(1, N), rep(0, T + 1 + N))
1946
        )
1947
        Rglpk.dir <- c(">=", "<=")
1948
        Rglpk.rhs <- c(min_sum, max_sum)
1949
  
1950
        # box constraint
1951
        Rglpk.bound <- list(
1952
          lower = list(
1953
            ind = c(1:N, N + T + 1),
1954
            val = c(constraints$min, -Inf)
1955
          ),
1956
          upper = list(
1957
            ind = c(1:N, N + T + 1),
1958
            val = c(constraints$max, Inf)
1959
          )
1960
        )
1961
  
1962
        # group constraint
1963
        if (!is.null(constraints$groups)) {
1964
          # check group constraint validility
1965
          if (!all(c(
1966
            length(constraints$cLO),
1967
            length(constraints$cUP)
1968
          )
1969
          == length(constraints$groups))) {
1970
            stop("Please assign group constraint")
1971
          }
1972
  
1973
          # group weight constraint
1974
          for (i in 1:length(constraints$groups)) {
1975
            # temp index variable
1976
            temp <- rep(0, 2 * N + T + 1)
1977
            temp[constraints$groups[[i]]] <- 1
1978
  
1979
            # weight constraints
1980
            Rglpk.mat <- rbind(Rglpk.mat, temp, temp)
1981
            Rglpk.dir <- c(Rglpk.dir, ">=", "<=")
1982
            Rglpk.rhs <- c(
1983
              Rglpk.rhs,
1984
              constraints$cLO[i],
1985
              constraints$cUP[i]
1986
            )
1987
          }
1988
  
1989
          # group postion limitation
1990
          if (!is.null(constraints$group_pos)) {
1991
            # check group position limitation length
1992
            if (length(constraints$group_pos)
1993
            != length(constraints$groups)) {
1994
              stop("Please assign group constraint")
1995
            }
1996
            for (i in 1:length(constraints$groups)) {
1997
              # temp index variable
1998
              temp <- rep(0, 2 * N + T + 1)
1999
              temp[constraints$groups[[i]] + N + T + 1] <- 1
2000
  
2001
              # position limitation constraints
2002
              Rglpk.mat <- rbind(Rglpk.mat, temp)
2003
              Rglpk.dir <- c(Rglpk.dir, "==")
2004
              Rglpk.rhs <- c(
2005
                Rglpk.rhs,
2006
                constraints$group_pos[i]
2007
              )
2008
            }
2009
          }
2010
  
2011
          # remove temp variable
2012
          rm(temp)
2013
        }
2014
  
2015
        # return target constraint
2016
        if (!is.infinite(target)) {
2017
          Rglpk.mat <- rbind(
2018
            Rglpk.mat,
2019
            c(colMeans(R), rep(0, N + T + 1))
2020
          )
2021
          Rglpk.dir <- c(Rglpk.dir, ">=")
2022
          Rglpk.rhs <- c(Rglpk.rhs, target)
2023
        }
2024
  
2025
        # position limitation constraint
2026
        if (!is.null(constraints$max_pos)) {
2027
          Rglpk.mat <- rbind(
2028
            Rglpk.mat,
2029
            c(rep(0, N + T + 1), rep(1, N))
2030
          )
2031
          Rglpk.dir <- c(Rglpk.dir, "<=")
2032
          Rglpk.rhs <- c(Rglpk.rhs, constraints$max_pos)
2033
  
2034
          Rglpk.mat <- rbind(
2035
            Rglpk.mat,
2036
            cbind(
2037
              diag(1, N),
2038
              matrix(0, N, T + 1),
2039
              diag(-upper, N)
2040
            ),
2041
            cbind(
2042
              diag(1, N),
2043
              matrix(0, N, T + 1),
2044
              diag(-lower, N)
2045
            )
2046
          )
2047
          Rglpk.dir <- c(Rglpk.dir, rep("<=", N), rep(">=", N))
2048
          Rglpk.rhs <- c(Rglpk.rhs, rep(0, 2 * N))
2049
        }
2050
  
2051
        # scenario constraint
2052
        Rglpk.mat <- rbind(
2053
          Rglpk.mat,
2054
          cbind(
2055
            matrix(R, T),
2056
            diag(1, T),
2057
            rep(1, T),
2058
            matrix(0, T, N)
2059
          )
2060
        )
2061
        Rglpk.dir <- c(Rglpk.dir, rep(">=", T))
2062
        Rglpk.rhs <- c(Rglpk.rhs, rep(0, T))
2063
  
2064
        # result from solver
2065
        Rglpk.result <- try(Rglpk::Rglpk_solve_LP(
2066
          obj = Rglpk.obj,
2067
          mat = Rglpk.mat,
2068
          dir = Rglpk.dir,
2069
          rhs = Rglpk.rhs,
2070
          bound = Rglpk.bound,
2071
          types = c(rep("C", N + T + 1), rep("B", N)),
2072
          max = TRUE
2073
        ))
2074
  
2075
        # null result
2076
        if (inherits(Rglpk.result, "try-error")) Rglpk.result <- NULL
2077
        if (is.null(Rglpk.result)) {
2078
          message(paste("Optimizer was unable to find a solution for target"))
2079
          return(paste("Optimizer was unable to find a solution for target"))
2080
        }
2081
  
2082
        # prepare output
2083
        weights <- as.vector(Rglpk.result$solution[1:N])
2084
        weights <- normalize_weights(weights)
2085
        names(weights) <- colnames(R)
2086
        obj_vals <- constrained_objective(
2087
          w = weights, R = R, portfolio = portfolio,
2088
          trace = TRUE, env = dotargs
2089
        )$objective_measures
2090
        out <- list(
2091
          weights = weights,
2092
          objective_measures = obj_vals,
2093
          opt_values = obj_vals,
2094
          out = Rglpk.result$optimum,
2095
          call = call
2096
        )
2097
        if (isTRUE(trace)) {
2098
          out$Rglpkoutput <- Rglpk.result
2099
        }
2100
      }
2101
  
2102
      # min ratio
2103
      else if (risk & reward) {
2104
        # utility function coef: weight + loss + VaR + shrinkage + position
2105
        Rglpk.obj <- c(
2106
          rep(0, N), # weight
2107
          rep(-1 / alpha / T, T), # loss
2108
          -1, # VaR
2109
          0, # shrinkage
2110
          rep(0, N) # position index
2111
        )
2112
  
2113
        # leverage constraint
2114
        Rglpk.mat <- rbind(
2115
          c(
2116
            rep(1, N), rep(0, T + 1),
2117
            -min_sum, rep(0, N)
2118
          ),
2119
          c(
2120
            rep(1, N), rep(0, T + 1),
2121
            -max_sum, rep(0, N)
2122
          )
2123
        )
2124
        Rglpk.dir <- c(">=", "<=")
2125
        Rglpk.rhs <- c(0, 0)
2126
  
2127
        # box constraint
2128
        Rglpk.mat <- rbind(
2129
          Rglpk.mat,
2130
          cbind(diag(1, N), matrix(0, N, T + 1), -lower, diag(0, N)),
2131
          cbind(diag(1, N), matrix(0, N, T + 1), -upper, diag(0, N))
2132
        )
2133
        Rglpk.dir <- c(Rglpk.dir, rep(">=", N), rep("<=", N))
2134
        Rglpk.rhs <- c(Rglpk.rhs, rep(0, 2 * N))
2135
  
2136
        # box constraint
2137
        Rglpk.bound <- list(
2138
          lower = list(
2139
            ind = c(1:N, N + T + 1),
2140
            val = rep(-Inf, N + 1)
2141
          ),
2142
          upper = list(
2143
            ind = c(1:N, N + T + 1),
2144
            val = rep(Inf, N + 1)
2145
          )
2146
        )
2147
  
2148
        # group constraint
2149
        if (!is.null(constraints$groups)) {
2150
          # check group constraint validility
2151
          if (!all(c(
2152
            length(constraints$cLO),
2153
            length(constraints$cUP)
2154
          )
2155
          == length(constraints$groups))) {
2156
            stop("Please assign group constraint")
2157
          }
2158
  
2159
          # update constraints matrix, direction and right-hand vector
2160
          for (i in 1:length(constraints$groups)) {
2161
            # temp index variable
2162
            # low level
2163
            temp <- c(rep(0, N + T + 1), -constraints$cLO[i], rep(0, N))
2164
            temp[constraints$groups[[i]]] <- 1
2165
            Rglpk.mat <- rbind(Rglpk.mat, temp)
2166
            Rglpk.dir <- c(Rglpk.dir, ">=")
2167
            Rglpk.rhs <- c(Rglpk.rhs, 0)
2168
  
2169
            # up level
2170
            temp <- c(rep(0, N + T + 1), -constraints$cUP[i], rep(0, N))
2171
            temp[constraints$groups[[i]]] <- 1
2172
            Rglpk.mat <- rbind(Rglpk.mat, temp)
2173
            Rglpk.dir <- c(Rglpk.dir, "<=")
2174
            Rglpk.rhs <- c(Rglpk.rhs, 0)
2175
          }
2176
  
2177
          # update group position constraints
2178
          if (!is.null(constraints$groups_pos)) {
2179
            # check group constraint validility
2180
            if (length(constraints$groups_pos)
2181
            != length(constraints$groups)) {
2182
              stop("Please assign group constraint")
2183
            }
2184
            for (i in 1:length(constraints$groups)) {
2185
              # temp index variable
2186
              # low level
2187
              temp <- c(rep(0, N + T + 2), rep(1, N))
2188
              Rglpk.mat <- rbind(Rglpk.mat, temp)
2189
              Rglpk.dir <- c(Rglpk.dir, "==")
2190
              Rglpk.rhs <- c(Rglpk.rhs, constraints$groups_pos[i])
2191
            }
2192
          }
2193
          rm(temp)
2194
        }
2195
  
2196
        # return target constraint
2197
        if (!is.infinite(target)) {
2198
          Rglpk.mat <- rbind(
2199
            Rglpk.mat,
2200
            c(rep(0, N + T + 1)),
2201
            target,
2202
            rep(0, N)
2203
          )
2204
  
2205
          Rglpk.dir <- c(Rglpk.dir, "<=")
2206
          Rglpk.rhs <- c(Rglpk.rhs, 1)
2207
        }
2208
  
2209
        # shrinkage constraint
2210
        Rglpk.mat <- rbind(
2211
          Rglpk.mat,
2212
          c(colMeans(R), rep(0, T + N + 2))
2213
        )
2214
        Rglpk.dir <- c(Rglpk.dir, "==")
2215
        Rglpk.rhs <- c(Rglpk.rhs, 1)
2216
  
2217
        # position limitation constraint
2218
        if (!is.null(constraints$max_pos)) {
2219
          Rglpk.mat <- rbind(
2220
            Rglpk.mat,
2221
            c(rep(0, N + T + 2), rep(1, N))
2222
          )
2223
          Rglpk.dir <- c(Rglpk.dir, "<=")
2224
          Rglpk.rhs <- c(Rglpk.rhs, constraints$max_pos)
2225
  
2226
          Rglpk.mat <- rbind(
2227
            Rglpk.mat,
2228
            cbind(
2229
              diag(1, N),
2230
              matrix(0, N, T + 2),
2231
              diag(-upper - 9999, N)
2232
            ),
2233
            cbind(
2234
              diag(1, N),
2235
              matrix(0, N, T + 2),
2236
              diag(-lower + 9999, N)
2237
            )
2238
          )
2239
          Rglpk.dir <- c(Rglpk.dir, rep("<=", N), rep(">=", N))
2240
          Rglpk.rhs <- c(Rglpk.rhs, rep(0, N), rep(0, N))
2241
        }
2242
  
2243
        # scenario constraint
2244
        Rglpk.mat <- rbind(
2245
          Rglpk.mat,
2246
          cbind(
2247
            matrix(R, T),
2248
            diag(1, T),
2249
            rep(1, T),
2250
            matrix(0, T, N + 1)
2251
          )
2252
        )
2253
        Rglpk.dir <- c(Rglpk.dir, rep(">=", T))
2254
        Rglpk.rhs <- c(Rglpk.rhs, rep(0, T))
2255
  
2256
        # result from solver
2257
        Rglpk.result <- try(Rglpk::Rglpk_solve_LP(
2258
          obj = Rglpk.obj,
2259
          mat = Rglpk.mat,
2260
          dir = Rglpk.dir,
2261
          rhs = Rglpk.rhs,
2262
          bound = Rglpk.bound,
2263
          types = c(rep("C", N + T + 2), rep("B", N)),
2264
          max = TRUE
2265
        ))
2266
  
2267
        # null result
2268
        if (inherits(Rglpk.result, "try-error")) Rglpk.result <- NULL
2269
        if (is.null(Rglpk.result)) {
2270
          message(paste("Optimizer was unable to find a solution for target"))
2271
          return(paste("Optimizer was unable to find a solution for target"))
2272
        }
2273
  
2274
        # prepare output
2275
        weights <- as.vector(Rglpk.result$solution[1:N] / Rglpk.result$solution[N + T + 2])
2276
        weights <- normalize_weights(weights)
2277
        names(weights) <- colnames(R)
2278
        obj_vals <- constrained_objective(
2279
          w = weights, R = R, portfolio = portfolio,
2280
          trace = TRUE, env = dotargs
2281
        )$objective_measures
2282
        out <- list(
2283
          weights = weights,
2284
          objective_measures = obj_vals,
2285
          opt_values = obj_vals,
2286
          out = Rglpk.result$optimum,
2287
          call = call
2288
        )
2289
        if (isTRUE(trace)) {
2290
          out$Rglpkoutput <- Rglpk.result
2291
        }
2292
      }
2293
    }
2294
  } ## end case for Rglpk
2295
  
2296
  ## case if method = osqp -- Operator Splitting Quadratic Program
2297
  if (optimize_method == "osqp") {
2298
    # search for required package
2299
    stopifnot("package:osqp" %in% search() ||
2300
      requireNamespace("osqp", quietly = TRUE))
2301
  
2302
    # osqp solver can only solve quadratic programming problems
2303
    valid_objnames <- c(
2304
      "mean", "sigma", "StdDev", "var", "volatility",
2305
      "sd"
2306
    )
2307
  
2308
    # default risk and reward parameter
2309
    reward <- FALSE
2310
    risk <- FALSE
2311
    target <- -Inf
2312
  
2313
    # search invalid objectives
2314
    for (objective in portfolio$objectives) {
2315
      if (objective$enabled) {
2316
        if (!(objective$name %in% valid_objnames)) {
2317
          stop("osqp only solves mean and sd type business objectives, 
2318
                 choose a different optimize method.")
2319
        }
2320
  
2321
        # return target
2322
        target <- ifelse(!is.null(objective$target), objective$target, target)
2323
  
2324
        # optimization objective function
2325
        reward <- ifelse(objective$name == "mean", TRUE, reward)
2326
        risk <- ifelse(objective$name %in% valid_objnames[2:6], TRUE, risk)
2327
      }
2328
    }
2329
  
2330
    # trim target return
2331
    target <- max(target, constraints$return_target, na.rm = TRUE)
2332
  
2333
    # get leverage paramters from constraints
2334
    max_sum <- constraints$max_sum
2335
    min_sum <- constraints$min_sum
2336
  
2337
    # transfer inf to big number
2338
    max_sum <- ifelse(is.infinite(max_sum), 9999.0, max_sum)
2339
    min_sum <- ifelse(is.infinite(min_sum), -9999.0, min_sum)
2340
  
2341
    # get upper and lower limitation
2342
    upper <- constraints$max
2343
    lower <- constraints$min
2344
  
2345
    upper[which(is.infinite(upper))] <- 9999.0
2346
    lower[which(is.infinite(lower))] <- 9999.0
2347
  
2348
    # issue message if min_sum and max_sum are restrictive
2349
    if ((max_sum - min_sum) < 0.02) {
2350
      message("Leverage constraint min_sum and max_sum are restrictive, 
2351
                consider relaxing. e.g. 'full_investment' constraint 
2352
                should be min_sum=0.99 and max_sum=1.01")
2353
    }
2354
  
2355
    # implement for different case
2356
    if (xor(risk, reward)) {
2357
      # set P matrix and q vector
2358
      if (reward) {
2359
        osqp.P <- diag(0, N)
2360
        osqp.q <- -colMeans(R)
2361
      }
2362
      else {
2363
        osqp.P <- cov(R)
2364
        osqp.q <- rep(0, N)
2365
      }
2366
  
2367
      # leverage constraint
2368
      osqp.A <- rep(1, N)
2369
      osqp.l <- min_sum
2370
      osqp.u <- max_sum
2371
      osqp.rnames <- c("sum")
2372
  
2373
      # box constraint
2374
      osqp.A <- rbind(osqp.A, diag(1, N))
2375
      osqp.l <- c(osqp.l, lower)
2376
      osqp.u <- c(osqp.u, upper)
2377
      osqp.rnames <- c(
2378
        osqp.rnames,
2379
        rep("box", N)
2380
      )
2381
  
2382
      # group constraint
2383
      if (!is.null(constraints$groups)) {
2384
        # check group constraint validility
2385
        if (!all(c(
2386
          length(constraints$cLO),
2387
          length(constraints$cUP)
2388
        )
2389
        == length(constraints$groups))) {
2390
          stop("Please assign group constraint")
2391
        }
2392
  
2393
        # update constraints matrix, direction and right-hand vector
2394
        for (i in 1:length(constraints$groups)) {
2395
          # temp index variable
2396
          temp <- rep(0, N)
2397
          temp[constraints$groups[[i]]] <- 1
2398
  
2399
          osqp.A <- rbind(osqp.A, temp)
2400
          osqp.l <- c(osqp.l, constraints$cLO[i])
2401
          osqp.u <- c(osqp.u, constraints$cUP[i])
2402
        }
2403
        rm(temp)
2404
        osqp.rnames <- c(
2405
          osqp.rnames,
2406
          rep(
2407
            "group",
2408
            length(constraints$groups)
2409
          )
2410
        )
2411
      }
2412
  
2413
      # return target constraint
2414
      if (!is.infinite(target)) {
2415
        osqp.A <- rbind(osqp.A, colMeans(R))
2416
        osqp.l <- c(osqp.l, target)
2417
        osqp.u <- c(osqp.u, Inf)
2418
        osqp.rnames <- c(osqp.rnames, "target")
2419
      }
2420
  
2421
      rownames(osqp.A) <-
2422
        names(osqp.u) <-
2423
        names(osqp.l) <- osqp.rnames
2424
      colnames(osqp.A) <- colnames(R)
2425
  
2426
      # result from solver
2427
      osqp.result <- try(osqp::solve_osqp(
2428
        P = osqp.P,
2429
        q = osqp.q,
2430
        A = osqp.A,
2431
        l = osqp.l,
2432
        u = osqp.u,
2433
        pars = osqp::osqpSettings(verbose = FALSE)
2434
      ))
2435
  
2436
      # null result
2437
      if (inherits(osqp.result, "try-error")) osqp.result <- NULL
2438
      if (is.null(osqp.result)) {
2439
        message(paste("Optimizer was unable to find a solution for target"))
2440
        return(paste("Optimizer was unable to find a solution for target"))
2441
      }
2442
  
2443
      # prepare output
2444
      weights <- as.vector(osqp.result$x)
2445
      weights <- normalize_weights(weights)
2446
      names(weights) <- colnames(R)
2447
      obj_vals <- constrained_objective(
2448
        w = weights, R = R, portfolio = portfolio,
2449
        trace = TRUE, env = dotargs
2450
      )$objective_measures
2451
      out <- list(
2452
        weights = weights,
2453
        objective_measures = obj_vals,
2454
        opt_values = obj_vals,
2455
        out = osqp.result$info$obj_val,
2456
        call = call
2457
      )
2458
      if (isTRUE(trace)) {
2459
        out$osqpoutput <- osqp.result
2460
      }
2461
    }
2462
    else {
2463
      osqp.P <- cbind(rbind(cov(R), rep(0, N)), rep(0, N + 1))
2464
      osqp.q <- rep(0, N + 1)
2465
      osqp.rnames <- c()
2466
  
2467
      # leverage constraint
2468
      osqp.A <- c(rep(1, N), -min_sum)
2469
      osqp.A <- rbind(osqp.A, c(rep(-1, N), max_sum))
2470
      osqp.l <- c(0, 0)
2471
      osqp.u <- c(Inf, Inf)
2472
      osqp.rnames <- c("min.sum", "max.sum")
2473
  
2474
      # box constraint
2475
      osqp.A <- rbind(
2476
        osqp.A,
2477
        cbind(diag(1, N), -lower)
2478
      )
2479
      osqp.A <- rbind(
2480
        osqp.A,
2481
        cbind(diag(-1, N), upper)
2482
      )
2483
      osqp.l <- c(osqp.l, rep(0, 2 * N))
2484
      osqp.u <- c(osqp.u, rep(Inf, 2 * N))
2485
      osqp.rnames <- c(osqp.rnames, rep("Box", 2 * N))
2486
  
2487
      # group constraint
2488
      if (!is.null(constraints$groups)) {
2489
        # check group constraint validility
2490
        if (!all(c(
2491
          length(constraints$cLO),
2492
          length(constraints$cUP)
2493
        )
2494
        == length(constraints$groups))) {
2495
          stop("Please assign group constraint")
2496
        }
2497
  
2498
        # update constraints matrix, direction and right-hand vector
2499
        for (i in 1:length(constraints$groups)) {
2500
          # temp index variable
2501
          temp <- rep(0, N)
2502
          temp[constraints$groups[[i]]] <- 1
2503
  
2504
          osqp.A <- rbind(
2505
            osqp.A,
2506
            c(temp, -constraints$cLO[i])
2507
          )
2508
          osqp.A <- rbind(
2509
            osqp.A,
2510
            c(-temp, constraints$cUP[i])
2511
          )
2512
          osqp.l <- c(osqp.l, 0, 0)
2513
          osqp.u <- c(osqp.u, Inf, Inf)
2514
        }
2515
  
2516
        osqp.rnames <- c(osqp.rnames, rep(
2517
          "Group",
2518
          2 * length(constraints$groups)
2519
        ))
2520
        rm(temp)
2521
      }
2522
  
2523
      # return target constraint
2524
      if (!is.infinite(target)) {
2525
        osqp.A <- rbind(osqp.A, c(rep(0, N), target))
2526
        osqp.l <- c(osqp.l, -Inf)
2527
        osqp.u <- c(osqp.u, 1)
2528
        osqp.rnames <- c(osqp.rnames, "target")
2529
      }
2530
  
2531
      # shrinkage constraint
2532
      osqp.A <- rbind(osqp.A, c(colMeans(R), 0))
2533
      osqp.l <- c(osqp.l, 1)
2534
      osqp.u <- c(osqp.u, 1)
2535
      osqp.A <- rbind(osqp.A, c(rep(0, N), 1))
2536
      osqp.l <- c(osqp.l, 0)
2537
      osqp.u <- c(osqp.u, Inf)
2538
      osqp.rnames <- c(osqp.rnames, rep("Shrinkage", 2))
2539
  
2540
      rownames(osqp.A) <-
2541
        names(osqp.u) <-
2542
        names(osqp.l) <- osqp.rnames
2543
      colnames(osqp.A) <- c(colnames(R), "Shrinkage")
2544
  
2545
      # result from solver
2546
      osqp.result <- try(osqp::solve_osqp(
2547
        P = osqp.P,
2548
        q = osqp.q,
2549
        A = osqp.A,
2550
        l = osqp.l,
2551
        u = osqp.u,
2552
        pars = osqp::osqpSettings(verbose = FALSE)
2553
      ))
2554
  
2555
      # null result
2556
      if (inherits(osqp.result, "try-error")) osqp.result <- NULL
2557
      if (is.null(osqp.result)) {
2558
        message(paste("Optimizer was unable to find a solution for target"))
2559
        return(paste("Optimizer was unable to find a solution for target"))
2560
      }
2561
  
2562
      # prepare output
2563
      weights <- as.vector(osqp.result$x[1:N] / osqp.result$x[N + 1])
2564
      weights <- normalize_weights(weights)
2565
      names(weights) <- colnames(R)
2566
      obj_vals <- constrained_objective(
2567
        w = weights, R = R, portfolio = portfolio,
2568
        trace = TRUE, env = dotargs
2569
      )$objective_measures
2570
      out <- list(
2571
        weights = weights,
2572
        objective_measures = obj_vals,
2573
        opt_values = obj_vals,
2574
        out = osqp.result$info$obj_val,
2575
        call = call
2576
      )
2577
      if (isTRUE(trace)) {
2578
        out$osqpoutput <- osqp.result
2579
      }
2580
    }
2581
  } ## end case for osqp
2582
  
2583
  ## case if method = mco -- Multi-criteria Optimization
2584
  if (optimize_method == "mco") {
2585
    # utility function
2586
    mco.fn <- function(w) {
2587
      # default return index
2588
      mco.return <- FALSE
2589
  
2590
      # default risk index
2591
      mco.risk <- FALSE
2592
  
2593
      # default lambda
2594
      mco.lambda <- FALSE
2595
  
2596
      # default alpha
2597
      mco.alpha <- 0.05
2598
  
2599
      # search reward and risk in active objectives
2600
      for (objective in portfolio$objectives) {
2601
        if (objective$enabled) {
2602
          # grab reward
2603
          mco.return <- ifelse(objective$name == "mean",
2604
            t(w) %*% colMeans(R),
2605
            mco.return
2606
          )
2607
  
2608
          # grab risk
2609
          mco.risk <- switch(
2610
            objective$name,
2611
            "ES" = 1,
2612
            "CVaR" = 1,
2613
            "TVaR" = 1,
2614
            "sigma" = 2,
2615
            "volitility" = 2,
2616
            "StdDev" = 2,
2617
            "sd" = 2
2618
          )
2619
  
2620
          # check risk
2621
          mco.risk <- ifelse(is.null(mco.risk), FALSE, mco.risk)
2622
  
2623
          # grab alpha
2624
          mco.alpha <- switch(
2625
            is.null(objective$arguments[["p"]]) + 1,
2626
            objective$arguments[["p"]],
2627
            mco.alpha
2628
          )
2629
  
2630
          # check alpha
2631
          mco.alpha <- ifelse(mco.alpha > 0.5, 1 - mco.alpha, mco.alpha)
2632
  
2633
          # grab lambda
2634
          mco.lambda <- ifelse(is.null(objective$risk_aversion),
2635
            mco.lambda,
2636
            objective$risk_aversion
2637
          )
2638
        }
2639
      }
2640
  
2641
      # temp portfolio
2642
      mco.portfolio <- R %*% w
2643
  
2644
      # risk calculation
2645
      if (mco.risk == 1) {
2646
        mco.VaR <- quantile(mco.portfolio, mco.alpha)
2647
        mco.output.risk <- -mean(mco.portfolio[which(mco.portfolio <= mco.VaR)])
2648
      }
2649
      if (mco.risk == 2) {
2650
        mco.output.risk <- sd(mco.portfolio)
2651
        if (mco.lambda) {
2652
          mco.output.risk <- mco.output.risk - mco.lambda * mean(mco.portfolio)
2653
        }
2654
      }
2655
  
2656
      # negative return portfolio
2657
      if (mean(mco.portfolio) < 0) {
2658
        return(Inf)
2659
      }
2660
  
2661
      # output
2662
      if (mco.return) {
2663
        if (mco.risk) {
2664
          return(mco.output.risk / mean(mco.portfolio))
2665
        }
2666
        return(-mean(mco.portfolio))
2667
      }
2668
      return(mco.output.risk)
2669
    }
2670
  
2671
    # input dimension
2672
    mco.idim <- N
2673
  
2674
    # output dimension
2675
    mco.odim <- 1
2676
  
2677
    # constraints
2678
    mco.constraints <- function(w) {
2679
      result <- 0
2680
  
2681
      # leverage constraint
2682
      result <- result - 999 * max(sum(w) > constraints$max_sum, 0)
2683
      result <- result - 999 * max(constraints$min_sum > sum(w), 0)
2684
  
2685
      # position limitation
2686
      result <- result -
2687
        999 * max(sum(w != 0) > constraints$max_pos, 0, na.rm = TRUE)
2688
  
2689
      # return target
2690
      result <- result -
2691
        999 * max(t(w) %*% colMeans(R) < constraints$return_target,
2692
          0,
2693
          na.rm = TRUE
2694
        )
2695
  
2696
      # diversity constraint
2697
      result <- result -
2698
        999 * max(constraints$div_target > (1 - sum(w^2)),
2699
          0,
2700
          na.rm = TRUE
2701
        )
2702
  
2703
      # group constraint
2704
      if (!is.null(constraints$groups)) {
2705
        for (i in 1:length(constraints$groups)) {
2706
          result <- result -
2707
            999 * max(sum(w[constraints$groups[[i]]]) > constraints$cUP[i],
2708
              0,
2709
              na.rm = TRUE
2710
            )
2711
          result <- result -
2712
            999 * max(constraints$cLO[i] > sum(w[constraints$groups[[i]]]),
2713
              0,
2714
              na.rm = TRUE
2715
            )
2716
        }
2717
        if (!is.null(constraints$group_pos)) {
2718
          for (i in 1:length(constraints$groups)) {
2719
            temp <- w[constraints$groups[[i]]]
2720
            result <- result -
2721
              999 * max(sum(temp != 0) > constraints$group_pos[i],
2722
                0,
2723
                na.rm = TRUE
2724
              )
2725
          }
2726
        }
2727
      }
2728
  
2729
      return(result)
2730
    }
2731
  
2732
    # check group constraint
2733
    if (!is.null(constraints$groups)) {
2734
      # check group constraint validility
2735
      if (!all(c(
2736
        length(constraints$cLO),
2737
        length(constraints$cUP)
2738
      )
2739
      == length(constraints$groups))) {
2740
        stop("Please assign group constraint")
2741
      }
2742
      if (!is.null(constraints$group_pos)) {
2743
        if (length(constraints$group_pos)
2744
        != length(constraints$groups)) {
2745
          stop("Please assign group constraint")
2746
        }
2747
      }
2748
    }
2749
  
2750
    # result from solver
2751
    mco.result <- try(mco::nsga2(
2752
      fn = mco.fn,
2753
      idim = mco.idim,
2754
      odim = mco.odim,
2755
      constraints = mco.constraints,
2756
      cdim = 1,
2757
      lower.bounds = constraints$min,
2758
      upper.bounds = constraints$max
2759
    ))
2760
  
2761
    # null result
2762
    if (inherits(mco.result, "try-error")) mco.result <- NULL
2763
    if (is.null(mco.result)) {
2764
      message(paste("Optimizer was unable to find a solution for target"))
2765
      return(paste("Optimizer was unable to find a solution for target"))
2766
    }
2767
  
2768
    # prepare output
2769
    weights <- as.vector(mco.result$par[1, ])
2770
    weights <- normalize_weights(weights)
2771
    names(weights) <- colnames(R)
2772
    obj_vals <- constrained_objective(
2773
      w = weights, R = R, portfolio = portfolio,
2774
      trace = TRUE, env = dotargs
2775
    )$objective_measures
2776
    out <- list(
2777
      weights = weights,
2778
      objective_measures = obj_vals,
2779
      opt_values = obj_vals,
2780
      out = mco.result$value[1, 1],
2781
      call = call
2782
    )
2783
    if (isTRUE(trace)) {
2784
      out$mcooutput <- mco.result
2785
    }
2786
  } ## end case for mco
2787
  
2788
  ## case if method = CVXR
2789
  cvxr_solvers <- c("CBC", "GLPK", "GLPK_MI", "OSQP", "CPLEX", "SCS", "ECOS", "GUROBI", "MOSEK")
2790
  if (optimize_method == "CVXR" || optimize_method == "cvxr" || optimize_method %in% cvxr_solvers) {
2791
    ## search for required package
2792
    stopifnot("package:CVXR" %in% search() || requireNamespace("CVXR", quietly = TRUE))
2793
    if(optimize_method == "CVXR"|| optimize_method == "cvxr") cvxr_default=TRUE else cvxr_default=FALSE
2794
    
2795
    ## variables
2796
    X <- as.matrix(R)
2797
    wts <- CVXR::Variable(N)
2798
    z <- CVXR::Variable(T)
2799
    zeta <- CVXR::Variable(1)
2800
    t <- CVXR::Variable(1)
2801
    
2802
    # objective type
2803
    target = -Inf
2804
    reward <- FALSE
2805
    risk <- FALSE
2806
    risk_ES <- FALSE
2807
    risk_CSM <- FALSE
2808
    risk_HHI <- FALSE
2809
    maxSR <- FALSE
2810
    maxSTARR <- FALSE
2811
    ESratio <- FALSE
2812
    CSMratio <- FALSE
2813
    alpha <- 0.05
2814
    lambda <- 1
2815
    
2816
    valid_objnames <- c("mean", "var", "sd", "StdDev", "ES", "CVaR", "ETL", "CSM", "HHI", "hhi")
2817
    for (objective in portfolio$objectives) {
2818
      if (objective$enabled) {
2819
        if (!(objective$name %in% valid_objnames)) {
2820
          stop("CVXR only solves mean, var/sd/StdDev and ETL/ES/CVaR/CSM/EQS type business objectives, 
2821
                 choose a different optimize_method.")
2822
        }
2823
        # alpha
2824
        alpha <- ifelse(!is.null(objective$arguments[["p"]]), objective$arguments[["p"]], alpha)
2825
        lambda <- ifelse(!is.null(objective$risk_aversion), objective$risk_aversion, lambda)
2826
        
2827
        # return target
2828
        target <- ifelse(!is.null(objective$target), objective$target, target)
2829
        
2830
        # optimization objective function
2831
        reward <- ifelse(objective$name == "mean", TRUE, reward)
2832
        risk <- ifelse(objective$name %in% valid_objnames[2:4], TRUE, risk)
2833
        risk_ES <- ifelse(objective$name %in% valid_objnames[5:7], TRUE, risk_ES)
2834
        risk_CSM <- ifelse(objective$name %in% valid_objnames[8:8], TRUE, risk_CSM)
2835
        if(objective$name %in% valid_objnames[9:11]){
2836
          risk_HHI <- TRUE
2837
          lambda_hhi <- objective$conc_aversion
2838
        }
2839
        arguments <- objective$arguments
2840
      }
2841
    }
2842
    if(alpha > 0.5) alpha <- (1 - alpha)
2843
    
2844
    mean_value <- mout$mu
2845
    sigma_value <- mout$sigma
2846
    
2847
    # ef will be called while generating efficient frontier
2848
    if(hasArg(ef)) ef=match.call(expand.dots=TRUE)$ef else ef=FALSE
2849
    if(ef){
2850
      reward=FALSE
2851
      mean_idx <- which(unlist(lapply(portfolio$objectives, function(x) x$name)) == "mean")
2852
      return_target <- portfolio$objectives[[mean_idx]]$target
2853
    } else return_target=NULL
2854
    
2855
    if(reward & !risk & !risk_ES & !risk_CSM & !risk_HHI){
2856
      # max return
2857
      obj <- -t(mean_value) %*% wts
2858
      constraints_cvxr = list()
2859
      tmpname = "mean"
2860
    } else if(!reward & risk & !risk_ES & !risk_CSM & !risk_HHI){
2861
      # min var/std
2862
      obj <- CVXR::quad_form(wts, sigma_value)
2863
      constraints_cvxr = list()
2864
      tmpname = "StdDev"
2865
    } else if(!reward & risk & !risk_ES & !risk_CSM & risk_HHI){
2866
      # min HHI
2867
      obj <- CVXR::quad_form(wts, sigma_value) + lambda_hhi * CVXR::cvxr_norm(wts, 2) **2
2868
      constraints_cvxr = list()
2869
      tmpname = "HHI"
2870
    } else if(reward & risk & !risk_ES & !risk_CSM & !risk_HHI){
2871
      # mean-var/sharpe ratio
2872
      if(hasArg(maxSR)) maxSR=match.call(expand.dots=TRUE)$maxSR
2873
      
2874
      if(!maxSR){
2875
        # min mean-variance
2876
        obj <- CVXR::quad_form(wts, sigma_value) - (t(mean_value) %*% wts) / lambda
2877
        constraints_cvxr = list()
2878
        tmpname = "optimal value"
2879
      } else {
2880
        # max sharpe ratio
2881
        obj <- CVXR::quad_form(wts, sigma_value)
2882
        constraints_cvxr = list(t(mean_value) %*% wts == 1, sum(wts) >= 0)
2883
        tmpname = "Sharpe Ratio"
2884
      }
2885
    } else if(risk_ES & !risk_CSM){
2886
      # ES objectives
2887
      if(reward){
2888
        if(hasArg(maxSTARR)) maxSTARR=match.call(expand.dots=TRUE)$maxSTARR else maxSTARR=TRUE
2889
        if(hasArg(ESratio)) maxSTARR=match.call(expand.dots=TRUE)$ESratio else maxSTARR=maxSTARR
2890
      }
2891
      if(maxSTARR){
2892
        # max ES ratio
2893
        obj <- zeta + (1/(T*alpha)) * sum(z)
2894
        constraints_cvxr = list(z >= 0, 
2895
                                z >= -X %*% wts - zeta, 
2896
                                t(mean_value) %*% wts == 1,
2897
                                sum(wts) >= 0)
2898
        tmpname = "ES ratio"
2899
      } else {
2900
        # min ES
2901
        obj <- zeta + (1/(T*alpha)) * sum(z)
2902
        constraints_cvxr = list(z >= 0, z >= -X %*% wts - zeta)
2903
        tmpname = "ES"
2904
      }
2905
    } else if(!risk_ES & risk_CSM){
2906
      # CSM objectives
2907
      if(reward){
2908
        if(hasArg(CSMratio)) CSMratio=match.call(expand.dots=TRUE)$CSMratio else CSMratio=TRUE
2909
      }
2910
      if(CSMratio){
2911
        # max CSM ratio
2912
        obj <- zeta + (1/(alpha* sqrt(T))) * t
2913
        constraints_cvxr = list(z >= 0, 
2914
                                z >= -X %*% wts - zeta, 
2915
                                t(mean_value) %*% wts == 1,
2916
                                sum(wts) >= 0,
2917
                                t >= CVXR::p_norm(z, p=2))
2918
        tmpname = "CSM ratio"
2919
      } else {
2920
        # min CSM
2921
        obj <- zeta + (1/(alpha* sqrt(T))) * t
2922
        constraints_cvxr = list(z >= 0, z >= -X %*% wts - zeta,
2923
                                t >= CVXR::p_norm(z, p=2))
2924
        tmpname = "CSM"
2925
      }
2926
    } else { 
2927
      # wrong objective
2928
      stop("Wrong multiple objectives. CVXR only solves mean, var, or simple ES/CSM type business objectives, please reorganize the objectives.")
2929
    }
2930
    
2931
    # weight scale for maximizing return per unit risk
2932
    if(!maxSR & !maxSTARR & !CSMratio) weight_scale=1 else weight_scale=sum(wts)
2933
    
2934
    # constraint type
2935
    ## weight sum constraint
2936
    ### Problem of maximizing return per unit risk doesn't need weight sum constraint, 
2937
    ### because weights could be scaled proportionally.
2938
    if(!maxSR & !maxSTARR & !CSMratio){
2939
      if(!is.null(constraints$max_sum) & !is.infinite(constraints$max_sum) & constraints$max_sum - constraints$min_sum <= 0.001){
2940
        constraints_cvxr = append(constraints_cvxr, sum(wts) == constraints$max_sum)
2941
      } else{
2942
        if(!is.null(constraints$max_sum)){
2943
          max_sum <- ifelse(is.infinite(constraints$max_sum), 9999.0, constraints$max_sum)
2944
          constraints_cvxr = append(constraints_cvxr, sum(wts) <= max_sum)
2945
        }
2946
        if(!is.null(constraints$min_sum)){
2947
          min_sum <- ifelse(is.infinite(constraints$min_sum), -9999.0, constraints$min_sum)
2948
          constraints_cvxr = append(constraints_cvxr, sum(wts) >= min_sum)
2949
        }
2950
      }
2951
    }
2952
    
2953
    ## box constraint
2954
    upper <- constraints$max
2955
    lower <- constraints$min
2956
    if(!all(is.infinite(upper))){
2957
      upper[which(is.infinite(upper))] <- 9999.0
2958
      constraints_cvxr = append(constraints_cvxr, wts <= upper * weight_scale)
2959
    }
2960
    if(!all(is.infinite(lower))){
2961
      lower[which(is.infinite(lower))] <- -9999.0
2962
      constraints_cvxr = append(constraints_cvxr, wts >= lower * weight_scale)
2963
    }
2964
    
2965
    ## group constraint
2966
    i = 1
2967
    for(g in constraints$groups){
2968
      constraints_cvxr = append(constraints_cvxr, sum(wts[g]) >= constraints$cLO[i] * weight_scale)
2969
      constraints_cvxr = append(constraints_cvxr, sum(wts[g]) <= constraints$cUP[i] * weight_scale)
2970
      i = i + 1
2971
    }
2972
    
2973
    ## target return constraint
2974
    ### target return in constraint
2975
    if(!is.null(constraints$return_target)){
2976
      constraints_cvxr = append(constraints_cvxr, t(mean_value) %*% wts >= constraints$return_target * weight_scale)
2977
    }
2978
    ### target return in objective, which is called by ef functions.
2979
    if(!is.null(return_target)){
2980
      constraints_cvxr = append(constraints_cvxr, t(mean_value) %*% wts >= return_target)
2981
    }
2982
    
2983
    ## factor exposure constraint
2984
    if(!is.null(constraints$B)){
2985
      constraints_cvxr = append(constraints_cvxr, t(constraints$B) %*% wts <= constraints$upper)
2986
      constraints_cvxr = append(constraints_cvxr, t(constraints$B) %*% wts >= constraints$lower)
2987
    }
2988
    
2989
    ## turnover constraint
2990
    if(!is.null(constraints$turnover_target)){
2991
      # set weight initial
2992
      if(is.null(constraints$weight_initial)){
2993
        weight_initial <- rep(1/N, N)
2994
      } else {
2995
        weight_initial <- constraints$weight_initial
2996
      }
2997
      # penalty for minvar
2998
      if(tmpname == "StdDev"){
2999
        if(is.null(constraints$turnover_penalty)){
3000
          sigma_value_penalty = nearPD(sigma_value)$mat
3001
        } else {
3002
          sigma_value_penalty = sigma_value + diag(constraints$turnover_penalty, N)
3003
        }
3004
        obj <- CVXR::quad_form(wts - weight_initial, sigma_value_penalty) + 2 * t(wts - weight_initial) %*% sigma_value %*% weight_initial + t(weight_initial) %*% sigma_value %*% weight_initial
3005
      }
3006
      constraints_cvxr = append(constraints_cvxr, sum(abs(wts - weight_initial)) <= constraints$turnover_target)
3007
    }
3008
    
3009
    # problem
3010
    prob_cvxr <- CVXR::Problem(CVXR::Minimize(obj), constraints = constraints_cvxr)
3011
    
3012
    if(cvxr_default){
3013
      if((risk || maxSR) && !risk_HHI){
3014
        result_cvxr <- CVXR::solve(prob_cvxr, solver = "OSQP", ... = ...)
3015
      } else {
3016
        result_cvxr <- CVXR::solve(prob_cvxr, solver = "SCS", ... = ...)
3017
      }
3018
    } else {
3019
      result_cvxr <- CVXR::solve(prob_cvxr, solver = optimize_method, ... = ...)
3020
    }
3021
    
3022
    cvxr_wts <- result_cvxr$getValue(wts)
3023
    if(maxSR | maxSTARR | CSMratio) cvxr_wts <- cvxr_wts / sum(cvxr_wts)
3024
    cvxr_wts <- as.vector(cvxr_wts)
3025
    cvxr_wts <- normalize_weights(cvxr_wts)
3026
    names(cvxr_wts) <- colnames(R)
3027
    
3028
    obj_cvxr <- list()
3029
    if(reward & !risk & !risk_ES & !risk_CSM & !risk_HHI){ # max return
3030
      obj_cvxr[[tmpname]] <- -result_cvxr$value
3031
    } else if(!reward & risk & !risk_ES & !risk_CSM & !risk_HHI){ # min var/std
3032
      obj_cvxr[[tmpname]] <- sqrt(result_cvxr$value)
3033
    } else if(!reward & risk & !risk_ES & !risk_CSM & risk_HHI){ # min HHI
3034
      obj_cvxr[["StdDev"]] <- sqrt(t(cvxr_wts) %*% sigma_value %*% cvxr_wts)
3035
      obj_cvxr[[tmpname]] <- (result_cvxr$value - t(cvxr_wts) %*% sigma_value %*% cvxr_wts)/lambda_hhi
3036
    } else if(!maxSR & !maxSTARR & !CSMratio){ # mean-var/ES/CSM
3037
      obj_cvxr[[tmpname]] <- result_cvxr$value
3038
      if(reward & risk){
3039
        obj_cvxr[["mean"]] <- cvxr_wts %*% mean_value
3040
        obj_cvxr[["StdDev"]] <- sqrt(t(cvxr_wts) %*% sigma_value %*% cvxr_wts)
3041
      }
3042
    } else { # max return per unit risk
3043
      obj_cvxr[["mean"]] <- cvxr_wts %*% mean_value
3044
      if(maxSR){
3045
        obj_cvxr[["StdDev"]] <- sqrt(t(cvxr_wts) %*% sigma_value %*% cvxr_wts)
3046
        obj_cvxr[[tmpname]] <- obj_cvxr[["mean"]] / obj_cvxr[["StdDev"]]
3047
      } else if(maxSTARR){
3048
        obj_cvxr[["ES"]] <- result_cvxr$value / sum(result_cvxr$getValue(wts))
3049
        obj_cvxr[[tmpname]] <- obj_cvxr[["mean"]] / obj_cvxr[["ES"]]
3050
      } else if(CSMratio){
3051
        obj_cvxr[["CSM"]] <- result_cvxr$value / sum(result_cvxr$getValue(wts))
3052
        obj_cvxr[[tmpname]] <- obj_cvxr[["mean"]] / obj_cvxr[["CSM"]]
3053
      }
3054
    }
3055
    if(ef) obj_cvxr[["mean"]] <- cvxr_wts %*% mean_value
3056
    
3057
    out = list(weights = cvxr_wts,
3058
               objective_measures = obj_cvxr,
3059
               opt_values=obj_cvxr,
3060
               out = obj_cvxr[[tmpname]],
3061
               call = call,
3062
               solver = result_cvxr$solver,
3063
               moment_values = list(momentFun = moment_name,mu = mout$mu, sigma = mout$sigma))
3064
    
3065
    optimize_method = "CVXR"
3066
  }## end case for CVXR
3067
  
3068
  # Prepare for final object to return
3069
  end_t <- Sys.time()
3070
  # print(c("elapsed time:",round(end_t-start_t,2),":diff:",round(diff,2), ":stats: ", round(out$stats,4), ":targets:",out$targets))
3071
  if(message) message(c("elapsed time:", end_t-start_t))
3072
  out$portfolio <- portfolio
3073
  if(trace) out$R <- R
3074
  out$data_summary <- list(first=first(R), last=last(R))
3075
  out$elapsed_time <- end_t - start_t
3076
  out$end_t <- as.character(Sys.time())
3077
  # return a $regime element to indicate what regime portfolio used for
3078
  # optimize.portfolio. The regime information is used in extractStats and
3079
  # extractObjectiveMeasures
3080
  if(regime.switching){
3081
    out$regime <- regime.idx
3082
  }
3083
  class(out) <- c(paste("optimize.portfolio", optimize_method, sep='.'), "optimize.portfolio")
3084
  return(out)
3085
}
3086

3087
#' @rdname optimize.portfolio.rebalancing
3088
#' @name optimize.portfolio.rebalancing
3089
#' @export
3090
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, rolling_window=NULL)
3091
{
3092
    stopifnot("package:foreach" %in% search() || requireNamespace("foreach",quietly=TRUE))
3093
    start_t<-Sys.time()
3094
      
3095
    #store the call for later
3096
    call <- match.call()
3097
    if(optimize_method=="random"){
3098
        # call random_portfolios() with constraints and search_size to create matrix of portfolios
3099
        if(is.null(rp))
3100
            rp<-random_portfolios(rpconstraints=constraints,permutations=search_size)
3101
    } else {
3102
        rp=NULL
3103
    }
3104
    
3105
    # check for trailing_periods argument and set rolling_window equal to 
3106
    # trailing_periods for backwards compatibility
3107
    if(hasArg(trailing_periods)) {
3108
      trailing_periods=match.call(expand.dots=TRUE)$trailing_periods
3109
      rolling_window <- trailing_periods
3110
    }
3111
    
3112
    if(is.null(training_period)) {if(nrow(R)<36) training_period=nrow(R) else training_period=36}
3113
    if (is.null(rolling_window)){
3114
        # define the index endpoints of our periods
3115
        ep.i<-endpoints(R,on=rebalance_on)[which(endpoints(R, on = rebalance_on)>=training_period)]
3116
        # now apply optimize.portfolio to the periods, in parallel if available
3117
        ep <- ep.i[1]
3118
        out_list<-foreach::foreach(ep=iterators::iter(ep.i), .errorhandling='pass', .packages='PortfolioAnalytics') %dopar% {
3119
                    optimize.portfolio(R[1:ep,],constraints=constraints,optimize_method=optimize_method, search_size=search_size, trace=trace, rp=rp, parallel=FALSE, ...=...)
3120
                  }
3121
    } else {
3122
        # define the index endpoints of our periods
3123
        ep.i<-endpoints(R,on=rebalance_on)[which(endpoints(R, on = rebalance_on)>=training_period)]
3124
        # now apply optimize.portfolio to the periods, in parallel if available
3125
        out_list<-foreach::foreach(ep=iterators::iter(ep.i), .errorhandling='pass', .packages='PortfolioAnalytics') %dopar% {
3126
                    optimize.portfolio(R[(ifelse(ep-rolling_window>=1,ep-rolling_window,1)):ep,],constraints=constraints,optimize_method=optimize_method, search_size=search_size, trace=trace, rp=rp, parallel=FALSE, ...=...)
3127
                  }
3128
    }
3129
    names(out_list)<-index(R[ep.i])
3130
    
3131
    end_t<-Sys.time()
3132
    message(c("overall elapsed time:",end_t-start_t))
3133
    class(out_list)<-c("optimize.portfolio.rebalancing")
3134
    return(out_list)
3135
}
3136

3137
#' Portfolio Optimization with Rebalancing Periods
3138
#' 
3139
#' Portfolio optimization with support for rebalancing periods for 
3140
#' out-of-sample testing (i.e. backtesting)
3141
#' 
3142
#' @details
3143
#' Run portfolio optimization with periodic rebalancing at specified time periods. 
3144
#' Running the portfolio optimization with periodic rebalancing can help 
3145
#' refine the constraints and objectives by evaluating the out of sample
3146
#' performance of the portfolio based on historical data.
3147
#' 
3148
#' If both \code{training_period} and \code{rolling_window} are \code{NULL}, 
3149
#' then \code{training_period} is set to a default value of 36. 
3150
#' 
3151
#' If \code{training_period} is \code{NULL} and a \code{rolling_window} is 
3152
#' specified, then \code{training_period} is set to the value of 
3153
#' \code{rolling_window}.
3154
#' 
3155
#' The user should be aware of the following behavior when both 
3156
#' \code{training_period} and \code{rolling_window} are specified and have 
3157
#' different values
3158
#' \describe{
3159
#'   \item{\code{training_period < rolling_window}: }{For example, if you have 
3160
#'   \code{rolling_window=60}, \code{training_period=50}, and the periodicity 
3161
#'   of the data is the same as the rebalance frequency (i.e. monthly data with 
3162
#'   \code{rebalance_on="months")} then the returns data used in the optimization 
3163
#'   at each iteration are as follows:
3164
#'   \itemize{
3165
#'   \item 1: R[1:50,]
3166
#'   \item 2: R[1:51,]
3167
#'   \item ...
3168
#'   \item 11: R[1:60,]
3169
#'   \item 12: R[1:61,]
3170
#'   \item 13: R[2:62,]
3171
#'   \item ...
3172
#'   }
3173
#'   This results in a growing window for several optimizations initially while
3174
#'   the endpoint iterator (i.e. \code{[50, 51, ...]}) is less than the 
3175
#'   rolling window width.}
3176
#'   \item{\code{training_period > rolling_window}: }{The data used in the initial 
3177
#'   optimization is \code{R[(training_period - rolling_window):training_period,]}. 
3178
#'   This results in some of the data being "thrown away", i.e. periods 1 to 
3179
#'   \code{(training_period - rolling_window - 1)} are not used in the optimization.}
3180
#' }
3181
#' 
3182
#' This function is a essentially a wrapper around \code{optimize.portfolio} 
3183
#' and thus the discussion in the Details section of the 
3184
#' \code{\link{optimize.portfolio}} help file is valid here as well.
3185
#' 
3186
#' This function is massively parallel and requires the 'foreach' package. It
3187
#' is suggested to register a parallel backend.
3188
#' 
3189
#' @param R an xts, vector, matrix, data frame, timeSeries or zoo object of asset returns
3190
#' @param portfolio an object of type "portfolio" specifying the constraints 
3191
#' and objectives for the optimization
3192
#' @param constraints default NULL, a list of constraint objects
3193
#' @param objectives default NULL, a list of objective objects
3194
#' @param optimize_method one of "DEoptim", "random", "pso", "GenSA", or "ROI"
3195
#' @param search_size integer, how many portfolios to test, default 20,000
3196
#' @param trace TRUE/FALSE if TRUE will attempt to return additional 
3197
#' information on the path or portfolios searched
3198
#' @param \dots any other passthru parameters to \code{\link{optimize.portfolio}}
3199
#' @param rp a set of random portfolios passed into the function to prevent recalculation
3200
#' @param rebalance_on character string of period to rebalance on. See 
3201
#' \code{\link[xts]{endpoints}} for valid names.
3202
#' @param training_period an integer of the number of periods to use as 
3203
#' a training data in the front of the returns data
3204
#' @param rolling_window an integer of the width (i.e. number of periods)
3205
#' of the rolling window, the default of NULL will run the optimization 
3206
#' using the data from inception.
3207
#' @return a list containing the following elements
3208
#' \describe{
3209
#'   \item{\code{portfolio}:}{ The portfolio object.}
3210
#'   \item{\code{R}:}{ The asset returns.}
3211
#'   \item{\code{call}:}{ The function call.}
3212
#'   \item{\code{elapsed_time:}}{ The amount of time that elapses while the 
3213
#'   optimization is run.}
3214
#'   \item{\code{opt_rebalancing:}}{ A list of \code{optimize.portfolio} 
3215
#'   objects computed at each rebalancing period.}
3216
#' }
3217
#' @author Kris Boudt, Peter Carl, Brian G. Peterson
3218
#' @name optimize.portfolio.rebalancing
3219
#' @aliases optimize.portfolio.rebalancing optimize.portfolio.rebalancing_v1
3220
#' @seealso \code{\link{portfolio.spec}} \code{\link{optimize.portfolio}}
3221
#' @examples
3222
#' \dontrun{
3223
#' data(edhec)
3224
#' R <- edhec[,1:4]
3225
#' funds <- colnames(R)
3226
#' 
3227
#' portf <- portfolio.spec(funds)
3228
#' portf <- add.constraint(portf, type="full_investment")
3229
#' portf <- add.constraint(portf, type="long_only")
3230
#' portf <- add.objective(portf, type="risk", name="StdDev")
3231
#' 
3232
#' # Quarterly rebalancing with 5 year training period
3233
#' bt.opt1 <- optimize.portfolio.rebalancing(R, portf,
3234
#' optimize_method="ROI",
3235
#' rebalance_on="quarters",
3236
#' training_period=60)
3237
#' 
3238
#' # Monthly rebalancing with 5 year training period and 4 year rolling window
3239
#' bt.opt2 <- optimize.portfolio.rebalancing(R, portf,
3240
#' optimize_method="ROI",
3241
#' rebalance_on="months",
3242
#' training_period=60,
3243
#' rolling_window=48)
3244
#' }
3245
#' @export
3246
optimize.portfolio.rebalancing <- function(R, portfolio=NULL, constraints=NULL, objectives=NULL, optimize_method=c("DEoptim", "random", "ROI", "CVXR"), search_size=20000, trace=FALSE, ..., rp=NULL, rebalance_on=NULL, training_period=NULL, rolling_window=NULL)
3247
{
3248
  stopifnot("package:foreach" %in% search() || requireNamespace("foreach",quietly=TRUE))
3249
  stopifnot("package:iterators" %in% search() || requireNamespace("iterators",quietly=TRUE))
3250
  
3251
  # This is the case where the user has passed in a list of portfolio objects
3252
  # for the portfolio argument.
3253
  # Loop through the portfolio list and recursively call 
3254
  # optimize.portfolio.rebalancing. 
3255
  #Note that I return at the end of this block. I know it is not good practice
3256
  # to return before the end of a function, but I am not sure of another way
3257
  # to handle a list of portfolio objects with the recursive call to 
3258
  # optimize.portfolio. 
3259
  if(inherits(portfolio, "portfolio.list")){
3260
    n.portf <- length(portfolio)
3261
    opt.list <- vector("list", n.portf)
3262
    for(i in 1:length(opt.list)){
3263
      if(hasArg(message)) message=match.call(expand.dots=TRUE)$message else message=FALSE
3264
      if(message) cat("Starting optimization of portfolio ", i, "\n")
3265
      opt.list[[i]] <- optimize.portfolio.rebalancing(R=R, 
3266
                                                      portfolio=portfolio[[i]], 
3267
                                                      constraints=constraints, 
3268
                                                      objectives=objectives, 
3269
                                                      optimize_method=optimize_method, 
3270
                                                      search_size=search_size, 
3271
                                                      trace=trace, 
3272
                                                      ...=..., 
3273
                                                      rp=rp, 
3274
                                                      rebalance_on=rebalance_on, 
3275
                                                      training_period=training_period, 
3276
                                                      rolling_window=rolling_window)
3277
    }
3278
    out <- combine.optimizations(opt.list)
3279
    class(out) <- "opt.rebal.list"
3280
    ##### return here for portfolio.list because this is a recursive call
3281
    ##### for optimize.portfolio.rebalancing
3282
    return(out)
3283
  }
3284
  
3285
  # This is the case where the user has passed in a mult.portfolio.spec
3286
  # object for multiple layer portfolio optimization.
3287
  if(inherits(portfolio, "mult.portfolio.spec")){
3288
    # This function calls optimize.portfolio.rebalancing on each sub portfolio
3289
    # according to the given optimization parameters and returns an xts object
3290
    # representing the proxy returns of each sub portfolio.
3291
    R <- proxy.mult.portfolio(R=R, mult.portfolio=portfolio)
3292
    # The optimization is controlled by the constraints and objectives in the
3293
    # top level portfolio
3294
    portfolio <- portfolio$top.portfolio
3295
  }
3296
  
3297
  # Store the call to return later
3298
  call <- match.call()
3299
  
3300
  start_t<-Sys.time()
3301
  
3302
  if (!is.null(portfolio) & !is.portfolio(portfolio)){
3303
    stop("you must pass in an object of class 'portfolio' to control the optimization")
3304
  }
3305
  
3306
  if(hasArg(message)) message=match.call(expand.dots=TRUE)$message else message=FALSE
3307
  
3308
  # check for trailing_periods argument and set rolling_window equal to 
3309
  # trailing_periods for backwards compatibility
3310
  if(hasArg(trailing_periods)) {
3311
    trailing_periods=match.call(expand.dots=TRUE)$trailing_periods
3312
    rolling_window <- trailing_periods
3313
  }
3314
  
3315
  
3316
  # Check for constraints and objectives passed in separately outside of the portfolio object
3317
  if(!is.null(constraints)){
3318
    if(inherits(constraints, "v1_constraint")){
3319
      if(is.null(portfolio)){
3320
        # If the user has not passed in a portfolio, we will create one for them
3321
        tmp_portf <- portfolio.spec(assets=constraints$assets)
3322
      }
3323
      message("constraint object passed in is a 'v1_constraint' object, updating to v2 specification")
3324
      portfolio <- update_constraint_v1tov2(portfolio=tmp_portf, v1_constraint=constraints)
3325
      # print.default(portfolio)
3326
    }
3327
    if(!inherits(constraints, "v1_constraint")){
3328
      # Insert the constraints into the portfolio object
3329
      portfolio <- insert_constraints(portfolio=portfolio, constraints=constraints)
3330
    }
3331
  }
3332
  if(!is.null(objectives)){
3333
    # Insert the objectives into the portfolio object
3334
    portfolio <- insert_objectives(portfolio=portfolio, objectives=objectives)
3335
  }
3336
  
3337
  #store the call for later
3338
  call <- match.call()
3339
  if(length(optimize_method) == 2) optimize_method <- optimize_method[2] else optimize_method <- optimize_method[1]
3340
  if(optimize_method=="random"){
3341
    # get any rp related arguments passed in through dots
3342
    if(hasArg(rp_method)) rp_method=match.call(expand.dots=TRUE)$rp_method else rp_method="sample"
3343
    if(hasArg(eliminate)) eliminate=match.call(expand.dots=TRUE)$eliminate else eliminate=TRUE
3344
    if(hasArg(fev)) fev=match.call(expand.dots=TRUE)$fev else fev=0:5
3345
    
3346
    # call random_portfolios() with constraints and search_size to create matrix of portfolios
3347
    if(is.null(rp))
3348
      if(inherits(portfolio, "regime.portfolios")){
3349
        rp <- rp.regime.portfolios(regime=portfolio, permutations=search_size, rp_method=rp_method, eliminate=eliminate, fev=fev)
3350
      } else {
3351
        rp <- random_portfolios(portfolio=portfolio, permutations=search_size, rp_method=rp_method, eliminate=eliminate, fev=fev)
3352
      }
3353
  } else {
3354
    rp = NULL
3355
  }
3356
  # set training_period equal to rolling_window if training_period is null
3357
  # and rolling_window is not null
3358
  if(is.null(training_period) & !is.null(rolling_window))
3359
    training_period <- rolling_window
3360
  
3361
  if(is.null(training_period)) {if(nrow(R)<36) training_period=nrow(R) else training_period=36}
3362
  
3363
  # turnover weight_initial is previous time point optimal weight
3364
  turnover_idx <- which(sapply(portfolio$constraints, function(x) x$type == "turnover"))
3365
  if(length(turnover_idx) > 0){
3366
    turnover_idx <- turnover_idx[1]
3367
    # original weight_initial
3368
    if(is.null(portfolio$constraints[[turnover_idx]]$weight_initial)){
3369
      portfolio$constraints[[turnover_idx]]$weight_initial <- rep(1/length(portfolio$assets), length(portfolio$assets))
3370
    }
3371
    
3372
    # rebalancing
3373
    if (is.null(rolling_window)){
3374
      # define the index endpoints of our periods
3375
      ep.i<-endpoints(R,on=rebalance_on)[which(endpoints(R, on = rebalance_on)>=training_period)]
3376
      # now apply optimize.portfolio to the periods, in parallel if available
3377
      ep <- ep.i[1]
3378
      out_list<-foreach::foreach(ep=iterators::iter(ep.i), .errorhandling='pass', .packages='PortfolioAnalytics') %dopar% {
3379
        opt <- optimize.portfolio(R[1:ep,], portfolio=portfolio, optimize_method=optimize_method, search_size=search_size, trace=trace, rp=rp, parallel=FALSE, ...=...)
3380
        portfolio$constraints[[turnover_idx]]$weight_initial <- opt$weights
3381
        opt
3382
      }
3383
    } else {
3384
      # define the index endpoints of our periods
3385
      ep.i<-endpoints(R,on=rebalance_on)[which(endpoints(R, on = rebalance_on)>=training_period)]
3386
      # now apply optimize.portfolio to the periods, in parallel if available
3387
      out_list<-foreach::foreach(ep=iterators::iter(ep.i), .errorhandling='pass', .packages='PortfolioAnalytics') %dopar% {
3388
        opt <- optimize.portfolio(R[(ifelse(ep-rolling_window>=1,ep-rolling_window,1)):ep,], portfolio=portfolio, optimize_method=optimize_method, search_size=search_size, trace=trace, rp=rp, parallel=FALSE, ...=...)
3389
        portfolio$constraints[[turnover_idx]]$weight_initial <- opt$weights
3390
        opt
3391
      }
3392
    }
3393
  } else {
3394
    if (is.null(rolling_window)){
3395
      # define the index endpoints of our periods
3396
      ep.i<-endpoints(R,on=rebalance_on)[which(endpoints(R, on = rebalance_on)>=training_period)]
3397
      # now apply optimize.portfolio to the periods, in parallel if available
3398
      ep <- ep.i[1]
3399
      out_list<-foreach::foreach(ep=iterators::iter(ep.i), .errorhandling='pass', .packages='PortfolioAnalytics') %dopar% {
3400
        optimize.portfolio(R[1:ep,], portfolio=portfolio, optimize_method=optimize_method, search_size=search_size, trace=trace, rp=rp, parallel=FALSE, ...=...)
3401
      }
3402
    } else {
3403
      # define the index endpoints of our periods
3404
      ep.i<-endpoints(R,on=rebalance_on)[which(endpoints(R, on = rebalance_on)>=training_period)]
3405
      # now apply optimize.portfolio to the periods, in parallel if available
3406
      out_list<-foreach::foreach(ep=iterators::iter(ep.i), .errorhandling='pass', .packages='PortfolioAnalytics') %dopar% {
3407
        optimize.portfolio(R[(ifelse(ep-rolling_window>=1,ep-rolling_window,1)):ep,], portfolio=portfolio, optimize_method=optimize_method, search_size=search_size, trace=trace, rp=rp, parallel=FALSE, ...=...)
3408
      }
3409
    }
3410
  }
3411
  
3412
  # out_list is a list where each element is an optimize.portfolio object
3413
  # at each rebalance date
3414
  names(out_list)<-index(R[ep.i])
3415
  
3416
  end_t <- Sys.time()
3417
  elapsed_time <- end_t - start_t
3418
  if(message) message(c("overall elapsed time:", end_t-start_t))
3419
  
3420
  # out object to return
3421
  out <- list()
3422
  out$portfolio <- portfolio
3423
  out$R <- R
3424
  out$call <- call
3425
  out$elapsed_time <- elapsed_time
3426
  out$opt_rebalancing <- out_list
3427
  
3428
  class(out) <- c("optimize.portfolio.rebalancing")
3429
  return(out)
3430
}
3431

3432
#' Execute multiple optimize.portfolio calls, presumably in parallel
3433
#' 
3434
#' This function will not speed up optimization!
3435
#' 
3436
#' This function exists to run multiple copies of optimize.portfolio, presumabley in parallel using foreach.
3437
#' 
3438
#' This is typically done to test your parameter settings, specifically 
3439
#' total population size, but also possibly to help tune your 
3440
#' convergence settings, number of generations, stopping criteria,
3441
#' etc.
3442
#' 
3443
#' If you want to use all the cores on your multi-core computer, use 
3444
#' the parallel version of the apppropriate optimization engine, not 
3445
#' this function.
3446
#' 
3447
#' @param R an xts, vector, matrix, data frame, timeSeries or zoo object of asset returns
3448
#' @param portfolio an object of type "portfolio" specifying the constraints and objectives for the optimization
3449
#' @param optimize_method one of "DEoptim", "random", "pso", "GenSA".
3450
#' @param search_size integer, how many portfolios to test, default 20,000
3451
#' @param trace TRUE/FALSE if TRUE will attempt to return additional information on the path or portfolios searched
3452
#' @param \dots any other passthru parameters
3453
#' @param rp matrix of random portfolio weights, default NULL, mostly for automated use by rebalancing optimization or repeated tests on same portfolios
3454
#' @param momentFUN the name of a function to call to set portfolio moments, default \code{\link{set.portfolio.moments_v2}}
3455
#' @param message TRUE/FALSE. The default is message=FALSE. Display messages if TRUE.
3456
#' @param nodes how many processes to run in the foreach loop, default 4
3457
#' 
3458
#' @return a list containing the optimal weights, some summary statistics, the function call, and optionally trace information 
3459
#' @author Kris Boudt, Peter Carl, Brian G. Peterson
3460
#' @export
3461
optimize.portfolio.parallel <- function(R,
3462
                                        portfolio,
3463
                                        optimize_method=c("DEoptim","random","ROI","pso","GenSA","CVXR"),
3464
                                        search_size=20000,
3465
                                        trace=FALSE, ...,
3466
                                        rp=NULL,
3467
                                        momentFUN='set.portfolio.moments',
3468
                                        message=FALSE,
3469
                                        nodes=4)
3470
{
3471
    stopifnot("package:foreach" %in% search() || requireNamespace("foreach",quietly=TRUE))
3472
    optimize_method=optimize_method[1]  
3473
    
3474
    start_t <- Sys.time()
3475
    
3476
    #store the call for later
3477
    call <- match.call()
3478
    
3479
    opt_out_list <- foreach::foreach(1:nodes, .errorhandling='pass', .packages='PortfolioAnalytics') %dopar% {
3480
      optimize.portfolio(R=R, portfolio=portfolio, 
3481
                         optimize_method=optimize_method, 
3482
                         search_size=search_size, trace=trace, 
3483
                         rp=rp, momentFUN=momentFUN, parallel=FALSE, 
3484
                         ...=...)
3485
    }
3486

3487
    end_t <- Sys.time()
3488
    elapsed_t <- end_t - start_t
3489
    if(message) message(c("overall elapsed time:", elapsed_t))
3490
    
3491
    out <- list()
3492
    out$optimizations <- opt_out_list
3493
    out$call <- call
3494
    out$elapsed_time <- elapsed_t
3495
    
3496
    class(out) <- c("optimize.portfolio.parallel")
3497
    return(out)
3498
}
3499

3500

3501
#TODO write function to compute an efficient frontier of optimal portfolios
3502

3503
###############################################################################
3504
# $Id$
3505
###############################################################################
3506

3507
Product

Resources

Company
