CoCalc -- optFUN.R

GitHub Repository: braverock/portfolioanalytics
Path: blob/master/R/optFUN.R
¹⁴³³ views
1

2
##### GMV and QU QP Function #####
3
#' GMV/QU QP Optimization
4
#' 
5
#' This function is called by optimize.portfolio to solve minimum variance or 
6
#' maximum quadratic utility problems
7
#' 
8
#' @param R xts object of asset returns
9
#' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
10
#' @param moments object of moments computed based on objective functions
11
#' @param lambda risk_aversion parameter
12
#' @param target target return value
13
#' @param lambda_hhi concentration aversion parameter
14
#' @param conc_groups list of vectors specifying the groups of the assets.
15
#' @param solver solver to use
16
#' @param control list of solver control parameters
17
#' @author Ross Bennett
18
gmv_opt <- function(R, constraints, moments, lambda, target, lambda_hhi, conc_groups, solver="quadprog", control=NULL){
19
  stopifnot("package:ROI" %in% search() || requireNamespace("ROI", quietly=TRUE))
20
  plugin <- paste0("ROI.plugin.", solver)
21
  stopifnot(paste0("package:", plugin) %in% search() || requireNamespace(plugin, quietly=TRUE))
22
  
23
  # Check for cleaned returns in moments
24
  if(!is.null(moments$cleanR)) R <- moments$cleanR
25
  
26
  # Number of assets
27
  N <- ncol(R)
28
  
29
  # Check for a target return constraint
30
  if(!is.na(target)) {
31
    # If var is the only objective specified, then moments$mean won't be calculated
32
    if(all(moments$mean==0)){
33
      tmp_means <- colMeans(R)
34
    } else {
35
      tmp_means <- moments$mean
36
    }
37
  } else {
38
    tmp_means <- rep(0, N)
39
    target <- 0
40
  }
41
  Amat <- tmp_means
42
  dir.vec <- "=="
43
  rhs.vec <- target
44
  meq <- 1
45
  
46
  # Set up initial A matrix for leverage constraints
47
  Amat <- rbind(Amat, rep(1, N), rep(-1, N))
48
  dir.vec <- c(dir.vec, ">=",">=")
49
  rhs.vec <- c(rhs.vec, constraints$min_sum, -constraints$max_sum)
50
  
51
  # Add min box constraints
52
  Amat <- rbind(Amat, diag(N))
53
  dir.vec <- c(dir.vec, rep(">=", N))
54
  rhs.vec <- c(rhs.vec, constraints$min)
55
  
56
  # Add max box constraints
57
  Amat <- rbind(Amat, -1*diag(N))
58
  dir.vec <- c(dir.vec, rep(">=", N))
59
  rhs.vec <- c(rhs.vec, -constraints$max)
60
  
61
  # Applying box constraints
62
  lb <- constraints$min
63
  ub <- constraints$max
64
  
65
  bnds <- ROI::V_bound(li=seq.int(1L, N), lb=as.numeric(lb),
66
                  ui=seq.int(1L, N), ub=as.numeric(ub))
67
  
68
  # Include group constraints
69
  if(try(!is.null(constraints$groups), silent=TRUE)){
70
    n.groups <- length(constraints$groups)
71
    Amat.group <- matrix(0, nrow=n.groups, ncol=N)
72
    for(i in 1:n.groups){
73
      Amat.group[i, constraints$groups[[i]]] <- 1
74
    }
75
    if(is.null(constraints$cLO)) cLO <- rep(-Inf, n.groups)
76
    if(is.null(constraints$cUP)) cUP <- rep(Inf, n.groups)
77
    Amat <- rbind(Amat, Amat.group, -Amat.group)
78
    dir.vec <- c(dir.vec, rep(">=", (n.groups + n.groups)))
79
    rhs.vec <- c(rhs.vec, constraints$cLO, -constraints$cUP)
80
  }
81
  
82
  # Add the factor exposures to Amat, dir.vec, and rhs.vec
83
  if(!is.null(constraints$B)){
84
    t.B <- t(constraints$B)
85
    Amat <- rbind(Amat, t.B, -t.B)
86
    dir.vec <- c(dir.vec, rep(">=", 2 * nrow(t.B)))
87
    rhs.vec <- c(rhs.vec, constraints$lower, -constraints$upper)
88
  }
89

90
  # quadprog cannot handle infinite values so replace Inf with .Machine$double.xmax
91
  # This is the strategy used in ROI
92
  # Amat[ is.infinite(Amat) & (Amat <= 0) ] <- -.Machine$double.xmax
93
  # Amat[ is.infinite(Amat) & (Amat >= 0) ] <-  .Machine$double.xmax
94
  # rhs.vec[is.infinite(rhs.vec) & (rhs.vec <= 0)] <- -.Machine$double.xmax
95
  # rhs.vec[is.infinite(rhs.vec) & (rhs.vec >= 0)] <- .Machine$double.xmax
96
  
97
  # Remove the rows of Amat and elements of dir.vec and rhs.vec where rhs.vec is Inf or -Inf
98
  Amat <- Amat[!is.infinite(rhs.vec), ]
99
  dir.vec <- dir.vec[!is.infinite(rhs.vec)]
100
  rhs.vec <- rhs.vec[!is.infinite(rhs.vec)]
101
  
102
  # Set up the quadratic objective
103
  if(!is.null(lambda_hhi)){
104
    if(length(lambda_hhi) == 1 & is.null(conc_groups)){
105
      ROI_objective <- ROI::Q_objective(Q=2*lambda*(moments$var + lambda_hhi * diag(N)), L=-moments$mean) # ROI
106
      #Dmat <- 2*lambda*(moments$var + lambda_hhi * diag(N)) # solve.QP
107
      #dvec <- moments$mean # solve.QP
108
    } else if(!is.null(conc_groups)){
109
      # construct the matrix with concentration aversion values by group
110
      hhi_mat <- matrix(0, nrow=N, ncol=N)
111
      vec <- 1:N
112
      for(i in 1:length(conc_groups)){
113
        tmpI <- diag(N)
114
        tmpvec <- conc_groups[[i]]
115
        zerovec <- setdiff(vec, tmpvec)
116
        for(j in 1:length(zerovec)){
117
          tmpI[zerovec[j], ] <- rep(0, N)
118
        }
119
        hhi_mat <- hhi_mat + lambda_hhi[i] * tmpI
120
      }
121
      ROI_objective <- ROI::Q_objective(Q=2*lambda*(moments$var + hhi_mat), L=-moments$mean) # ROI
122
      #Dmat <- 2 * lambda * (moments$var + hhi_mat) # solve.QP
123
      #dvec <- moments$mean # solve.QP
124
    }
125
  } else {
126
    ROI_objective <- ROI::Q_objective(Q=2*lambda*moments$var, L=-moments$mean) # ROI
127
    #Dmat <- 2 * lambda * moments$var # solve.QP
128
    #dvec <- moments$mean # solve.QP
129
  }
130
  # set up the optimization problem and solve
131
  opt.prob <- ROI::OP(objective=ROI_objective, 
132
                       constraints=ROI::L_constraint(L=Amat, dir=dir.vec, rhs=rhs.vec),
133
                 bounds=bnds)
134
  result <- try(ROI::ROI_solve(x=opt.prob, solver=solver, control=control), silent=TRUE)
135
  
136
  # result <- try(solve.QP(Dmat=Dmat, dvec=dvec, Amat=t(Amat), bvec=rhs.vec, meq=meq), silent=TRUE)
137
  if(inherits(x=result, "try-error")) stop(paste("No solution found:", result))
138
  
139
  weights <- result$solution[1:N]
140
  names(weights) <- colnames(R)
141
  out <- list()
142
  out$weights <- weights
143
  out$out <- result$objval
144
  obj_vals <- list()
145
  # Calculate the objective values here so that we can use the moments$mean
146
  # and moments$var that might be passed in by the user. This will avoid
147
  # the extra call to constrained_objective
148
  if(!all(moments$mean == 0)){
149
    port.mean <- as.numeric(sum(weights * moments$mean))
150
    names(port.mean) <- "mean"
151
    obj_vals[["mean"]] <- port.mean
152
    # faster and more efficient way to compute t(w) %*% Sigma %*% w
153
    port.sd <- sqrt(sum(crossprod(weights, moments$var) * weights))
154
    names(port.sd) <- "StdDev"
155
    obj_vals[["StdDev"]] <- port.sd
156
  } else {
157
    # faster and more efficient way to compute t(w) %*% Sigma %*% w
158
    port.sd <- sqrt(sum(crossprod(weights, moments$var) * weights))
159
    names(port.sd) <- "StdDev"
160
    obj_vals[["StdDev"]] <- port.sd
161
  }
162
  out$obj_vals <- obj_vals
163
  # out$out <- result$objval # ROI
164
  # out$call <- call # need to get the call outside of the function
165
  return(out)
166
}
167

168
##### Maximize Return LP Function #####
169
#' Maximum Return LP Optimization
170
#' 
171
#' This function is called by optimize.portfolio to solve maximum return
172
#' 
173
#' @param R xts object of asset returns
174
#' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
175
#' @param moments object of moments computed based on objective functions
176
#' @param target target return value
177
#' @param solver solver to use
178
#' @param control list of solver control parameters
179
#' @author Ross Bennett
180
maxret_opt <- function(R, moments, constraints, target, solver="glpk", control=NULL){
181
  stopifnot("package:ROI" %in% search() || requireNamespace("ROI",quietly = TRUE))
182
  plugin <- paste0("ROI.plugin.", solver)
183
  stopifnot(paste0("package:", plugin) %in% search() || requireNamespace(plugin, quietly=TRUE))
184
  
185
  # Check for cleaned returns in moments
186
  if(!is.null(moments$cleanR)) R <- moments$cleanR
187
  
188
  N <- ncol(R)
189
  # Applying box constraints
190
  # maxret_opt needs non infinite values for upper and lower bounds
191
  lb <- constraints$min
192
  ub <- constraints$max
193
  if(any(is.infinite(lb)) | any(is.infinite(ub))){
194
    warning("Inf or -Inf values detected in box constraints, maximum return 
195
            objectives must have finite box constraint values.")
196
    ub[is.infinite(ub)] <- max(abs(c(constraints$min_sum, constraints$max_sum)))
197
    lb[is.infinite(lb)] <- 0
198
  }
199
  bnds <- ROI::V_bound(li=seq.int(1L, N), lb=as.numeric(lb),
200
                  ui=seq.int(1L, N), ub=as.numeric(ub))
201
  
202
  # set up initial A matrix for leverage constraints
203
  Amat <- rbind(rep(1, N), rep(1, N))
204
  dir.vec <- c(">=","<=")
205
  rhs.vec <- c(constraints$min_sum, constraints$max_sum)
206
  
207
  # check for a target return
208
  if(!is.na(target)) {
209
    Amat <- rbind(Amat, moments$mean)
210
    dir.vec <- c(dir.vec, "==")
211
    rhs.vec <- c(rhs.vec, target)
212
  }
213
  
214
  # include group constraints
215
  if(try(!is.null(constraints$groups), silent=TRUE)){
216
    n.groups <- length(constraints$groups)
217
    Amat.group <- matrix(0, nrow=n.groups, ncol=N)
218
    for(i in 1:n.groups){
219
      Amat.group[i, constraints$groups[[i]]] <- 1
220
    }
221
    if(is.null(constraints$cLO)) cLO <- rep(-Inf, n.groups)
222
    if(is.null(constraints$cUP)) cUP <- rep(Inf, n.groups)
223
    Amat <- rbind(Amat, Amat.group, -Amat.group)
224
    dir.vec <- c(dir.vec, rep(">=", (n.groups + n.groups)))
225
    rhs.vec <- c(rhs.vec, constraints$cLO, -constraints$cUP)
226
  }
227
  
228
  # Add the factor exposures to Amat, dir.vec, and rhs.vec
229
  if(!is.null(constraints$B)){
230
    t.B <- t(constraints$B)
231
    Amat <- rbind(Amat, t.B, -t.B)
232
    dir.vec <- c(dir.vec, rep(">=", 2 * nrow(t.B)))
233
    rhs.vec <- c(rhs.vec, constraints$lower, -constraints$upper)
234
  }
235
  
236
  # set up the linear objective
237
  ROI_objective <- ROI::L_objective(L=-moments$mean)
238
  # objL <- -moments$mean
239
  
240
  # set up the optimization problem and solve
241
  opt.prob <- ROI::OP(objective=ROI_objective, 
242
                 constraints=ROI::L_constraint(L=Amat, dir=dir.vec, rhs=rhs.vec),
243
                 bounds=bnds)
244
  roi.result <- try(ROI::ROI_solve(x=opt.prob, solver=solver, control=control), silent=TRUE)
245
  if(inherits(roi.result, "try-error")) stop(paste("No solution found:", roi.result))
246
  
247
  # roi.result <- Rglpk_solve_LP(obj=objL, mat=Amat, dir=dir.vec, rhs=rhs.vec, bounds=bnds)
248
  
249
  # The Rglpk solvers status returns an an integer with status information
250
  # about the solution returned: 0 if the optimal solution was found, a 
251
  #non-zero value otherwise.
252
  if(roi.result$status$code != 0) {
253
    message(roi.result$status$msg$message)
254
    stop("No solution")
255
    return(NULL)
256
  }
257
  
258
  weights <- roi.result$solution[1:N]
259
  names(weights) <- colnames(R)
260
  out <- list()
261
  out$weights <- weights
262
  out$out <- roi.result$objval
263
  obj_vals <- list()
264
  # Calculate the objective values here so that we can use the moments$mean
265
  # that might be passed in by the user. This will avoid
266
  # the extra call to constrained_objective
267
  port.mean <- -roi.result$objval
268
  names(port.mean) <- "mean"
269
  obj_vals[["mean"]] <- port.mean
270
  out$obj_vals <- obj_vals
271
  # out$call <- call # need to get the call outside of the function
272
  return(out)
273
}
274

275
##### Maximize Return MILP Function #####
276
#' Maximum Return MILP Optimization
277
#' 
278
#' This function is called by optimize.portfolio to solve maximum return 
279
#' problems via mixed integer linear programming.
280
#' 
281
#' @param R xts object of asset returns
282
#' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
283
#' @param moments object of moments computed based on objective functions
284
#' @param target target return value
285
#' @param solver solver to use
286
#' @param control list of solver control parameters
287
#' @author Ross Bennett
288
maxret_milp_opt <- function(R, constraints, moments, target, solver="glpk", control=NULL){
289
  stopifnot("package:ROI" %in% search() || requireNamespace("ROI",quietly = TRUE))
290
  plugin <- paste0("ROI.plugin.", solver)
291
  stopifnot(paste0("package:", plugin) %in% search() || requireNamespace(plugin, quietly=TRUE))
292
  
293
  # Check for cleaned returns in moments
294
  if(!is.null(moments$cleanR)) R <- moments$cleanR
295
  
296
  # Number of assets
297
  N <- ncol(R)
298
  
299
  # Maximum number of positions (non-zero weights)
300
  max_pos <- constraints$max_pos
301
  min_pos <- 1
302
  
303
  # Upper and lower bounds on weights
304
  LB <- as.numeric(constraints$min)
305
  UB <- as.numeric(constraints$max)
306
  
307
  # Check for target return
308
  if(!is.na(target)){
309
    # We have a target
310
    targetcon <- rbind(c(moments$mean, rep(0, N)),
311
                       c(-moments$mean, rep(0, N)))
312
    targetdir <- c("<=", "==")
313
    targetrhs <- c(Inf, -target)
314
  } else {
315
    # No target specified, just maximize
316
    targetcon <- NULL
317
    targetdir <- NULL
318
    targetrhs <- NULL
319
  }
320
  
321
  # weight_sum constraint
322
  Amat <- rbind(c(rep(1, N), rep(0, N)),
323
                c(rep(1, N), rep(0, N)))
324
  
325
  # Target return constraint
326
  Amat <- rbind(Amat, targetcon)
327
  
328
  # Bounds and position limit constraints
329
  Amat <- rbind(Amat, cbind(-diag(N), diag(LB)))
330
  Amat <- rbind(Amat, cbind(diag(N), -diag(UB)))
331
  Amat <- rbind(Amat, c(rep(0, N), rep(-1, N)))
332
  Amat <- rbind(Amat, c(rep(0, N), rep(1, N)))
333
  
334
  dir <- c("<=", ">=", targetdir, rep("<=", 2*N), "<=", "<=")
335
  rhs <- c(1, 1, targetrhs, rep(0, 2*N), -min_pos, max_pos)
336
  
337
  # Include group constraints
338
  if(try(!is.null(constraints$groups), silent=TRUE)){
339
    n.groups <- length(constraints$groups)
340
    Amat.group <- matrix(0, nrow=n.groups, ncol=N)
341
    k <- 1
342
    l <- 0
343
    for(i in 1:n.groups){
344
      j <- constraints$groups[i] 
345
      Amat.group[i, k:(l+j)] <- 1
346
      k <- l + j + 1
347
      l <- k - 1
348
    }
349
    if(is.null(constraints$cLO)) cLO <- rep(-Inf, n.groups)
350
    if(is.null(constraints$cUP)) cUP <- rep(Inf, n.groups)
351
    zeros <- matrix(data=0, nrow=nrow(Amat.group), ncol=ncol(Amat.group))
352
    Amat <- rbind(Amat, cbind(Amat.group, zeros), cbind(-Amat.group, zeros))
353
    dir <- c(dir, rep(">=", (n.groups + n.groups)))
354
    rhs <- c(rhs, constraints$cLO, -constraints$cUP)
355
  }
356
  
357
  # Add the factor exposures to Amat, dir, and rhs
358
  if(!is.null(constraints$B)){
359
    t.B <- t(constraints$B)
360
    zeros <- matrix(data=0, nrow=nrow(t.B), ncol=ncol(t.B))
361
    Amat <- rbind(Amat, cbind(t.B, zeros), cbind(-t.B, zeros))
362
    dir <- c(dir, rep(">=", 2 * nrow(t.B)))
363
    rhs <- c(rhs, constraints$lower, -constraints$upper)
364
  }
365
  
366
  # Only seems to work if I do not specify bounds
367
  # bnds <- ROI::V_Bound(li=seq.int(1L, 2*m), lb=c(as.numeric(constraints$min), rep(0, m)),
368
  #                 ui=seq.int(1L, 2*m), ub=c(as.numeric(constraints$max), rep(1, m)))
369
  bnds <- NULL
370
  
371
  # Set up the types vector with continuous and binary variables
372
  types <- c(rep("C", N), rep("B", N))
373
  
374
  # Set up the linear objective to maximize mean return
375
  ROI_objective <- ROI::L_objective(L=c(-moments$mean, rep(0, N)))
376
  
377
  # Set up the optimization problem and solve
378
  opt.prob <- ROI::OP(objective=ROI_objective, 
379
                 constraints=ROI::L_constraint(L=Amat, dir=dir, rhs=rhs),
380
                 bounds=bnds, types=types)
381
  roi.result <- try(ROI::ROI_solve(x=opt.prob, solver=solver, control=control), silent=TRUE)
382
  if(inherits(roi.result, "try-error")) stop(paste("No solution found:", roi.result))
383
  
384
  # Weights
385
  weights <- roi.result$solution[1:N]
386
  names(weights) <- colnames(R)
387
  
388
  # The out object is returned
389
  out <- list()
390
  out$weights <- weights
391
  out$out <- roi.result$objval
392
  
393
  obj_vals <- list()
394
  port.mean <- -roi.result$objval
395
  names(port.mean) <- "mean"
396
  obj_vals[["mean"]] <- port.mean
397
  out$obj_vals <- obj_vals
398
  return(out)
399
}
400

401
##### Minimize ETL LP Function #####
402
#' Minimum ETL LP Optimization
403
#' 
404
#' This function is called by optimize.portfolio to solve minimum ETL problems.
405
#' 
406
#' @param R xts object of asset returns
407
#' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
408
#' @param moments object of moments computed based on objective functions
409
#' @param target target return value
410
#' @param alpha alpha value for ETL/ES/CVaR
411
#' @param solver solver to use
412
#' @param control list of solver control parameters
413
#' @author Ross Bennett
414
etl_opt <- function(R, constraints, moments, target, alpha, solver="glpk", control=NULL){
415
  stopifnot("package:ROI" %in% search() || requireNamespace("ROI",quietly = TRUE))
416
  plugin <- paste0("ROI.plugin.", solver)
417
  stopifnot(paste0("package:", plugin) %in% search() || requireNamespace(plugin, quietly=TRUE))
418
  
419
  # Check for cleaned returns in moments
420
  if(!is.null(moments$cleanR)) R <- moments$cleanR
421
  
422
  N <- ncol(R)
423
  T <- nrow(R)
424
  # Applying box constraints
425
  LB <- c(as.numeric(constraints$min), rep(0, T), -1)
426
  UB <- c(as.numeric(constraints$max), rep(Inf, T), 1)
427
  bnds <- ROI::V_bound(li=seq.int(1L, N+T+1), lb=LB,
428
                  ui=seq.int(1L, N+T+1), ub=UB)
429
  
430
  # Add this check if mean is not an objective and return is a constraints
431
  if(!is.na(target)){
432
    if(all(moments$mean == 0)){
433
      moments$mean <- colMeans(R)
434
    }
435
  } else {
436
    moments$mean <- rep(0, N)
437
    target <- 0
438
  }
439
  
440
  Amat <- cbind(rbind(1, 1, moments$mean, coredata(R)), rbind(0, 0, 0, cbind(diag(T), 1))) 
441
  dir.vec <- c(">=","<=",">=",rep(">=",T))
442
  rhs.vec <- c(constraints$min_sum, constraints$max_sum, target ,rep(0, T))
443
  
444
  if(try(!is.null(constraints$groups), silent=TRUE)){
445
    n.groups <- length(constraints$groups)
446
    Amat.group <- matrix(0, nrow=n.groups, ncol=N)
447
    for(i in 1:n.groups){
448
      Amat.group[i, constraints$groups[[i]]] <- 1
449
    }
450
    if(is.null(constraints$cLO)) cLO <- rep(-Inf, n.groups)
451
    if(is.null(constraints$cUP)) cUP <- rep(Inf, n.groups)
452
    zeros <- matrix(0, nrow=n.groups, ncol=(T+1))
453
    Amat <- rbind(Amat, cbind(Amat.group, zeros), cbind(-Amat.group, zeros))
454
    dir.vec <- c(dir.vec, rep(">=", (n.groups + n.groups)))
455
    rhs.vec <- c(rhs.vec, constraints$cLO, -constraints$cUP)
456
  }
457
  # Add the factor exposures to Amat, dir, and rhs
458
  if(!is.null(constraints$B)){
459
    t.B <- t(constraints$B)
460
    zeros <- matrix(data=0, nrow=nrow(t.B), ncol=(T+1))
461
    Amat <- rbind(Amat, cbind(t.B, zeros), cbind(-t.B, zeros))
462
    dir.vec <- c(dir.vec, rep(">=", 2 * nrow(t.B)))
463
    rhs.vec <- c(rhs.vec, constraints$lower, -constraints$upper)
464
  }
465
  
466
  ROI_objective <- ROI::L_objective(c(rep(0,N), rep(1/(alpha*T),T), 1))
467
  opt.prob <- ROI::OP(objective=ROI_objective, 
468
                       constraints=ROI::L_constraint(L=Amat, dir=dir.vec, rhs=rhs.vec),
469
                       bounds=bnds)
470
  roi.result <- try(ROI::ROI_solve(x=opt.prob, solver=solver, control=control), silent=TRUE)
471
  if(inherits(x=roi.result, "try-error")) stop(paste("No solution found:", roi.result))
472
  
473
  weights <- roi.result$solution[1:N]
474
  names(weights) <- colnames(R)
475
  out <- list()
476
  out$weights <- weights
477
  out$out <- roi.result$objval
478
  es_names <- c("ES", "ETL", "CVaR")
479
  es_idx <- which(es_names %in% names(moments))
480
  obj_vals <- list()
481
  # Calculate the objective values here so that we can use the moments$mean
482
  # and moments$var that might be passed in by the user. This will avoid
483
  # the extra call to constrained_objective
484
  if(!all(moments$mean == 0)){
485
    port.mean <- as.numeric(sum(weights * moments$mean))
486
    names(port.mean) <- "mean"
487
    obj_vals[["mean"]] <- port.mean
488
    port.es <- roi.result$objval
489
    names(port.es) <- es_names[es_idx]
490
    obj_vals[[es_names[es_idx]]] <- port.es
491
  } else {
492
    port.es <- roi.result$objval
493
    names(port.es) <- es_names[es_idx]
494
    obj_vals[[es_names[es_idx]]] <- port.es
495
  }
496
  out$obj_vals <- obj_vals
497
  return(out)
498
}
499

500
##### Minimize ETL MILP Function #####
501
#' Minimum ETL MILP Optimization
502
#' 
503
#' This function is called by optimize.portfolio to solve minimum ETL problems 
504
#' via mixed integer linear programming.
505
#' 
506
#' @param R xts object of asset returns
507
#' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
508
#' @param moments object of moments computed based on objective functions
509
#' @param target target return value
510
#' @param alpha alpha value for ETL/ES/CVaR
511
#' @param solver solver to use
512
#' @param control list of solver control parameters
513
#' @author Ross Bennett
514
etl_milp_opt <- function(R, constraints, moments, target, alpha, solver="glpk", control=NULL){
515
  stopifnot("package:ROI" %in% search() || requireNamespace("ROI",quietly = TRUE))
516
  plugin <- paste0("ROI.plugin.", solver)
517
  stopifnot(paste0("package:", plugin) %in% search() || requireNamespace(plugin, quietly=TRUE))
518
  
519
  # Check for cleaned returns in moments
520
  if(!is.null(moments$cleanR)) R <- moments$cleanR
521
  
522
  # Number of rows
523
  n <- nrow(R)
524
  
525
  # Number of columns
526
  m <- ncol(R)
527
  
528
  max_sum <- constraints$max_sum
529
  min_sum <- constraints$min_sum
530
  LB <- constraints$min
531
  UB <- constraints$max
532
  max_pos <- constraints$max_pos
533
  min_pos <- 1
534
  moments_mean <- as.numeric(moments$mean)
535
  
536
  # A benchmark can be specified in the parma package. 
537
  # Leave this in and set to 0 for now
538
  benchmark <- 0
539
  
540
  # Check for target return
541
  if(!is.na(target)){
542
    # We have a target
543
    targetcon <- c(moments_mean, rep(0, n+2))
544
    targetdir <- "=="
545
    targetrhs <- target
546
  } else {
547
    # No target specified, just maximize
548
    targetcon <- NULL
549
    targetdir <- NULL
550
    targetrhs <- NULL
551
  }
552
  
553
  # Set up initial A matrix
554
  tmpAmat <- cbind(-coredata(R),
555
                   matrix(-1, nrow=n, ncol=1), 
556
                   -diag(n),
557
                   matrix(benchmark, nrow=n, ncol=1))
558
  
559
  # Add leverage constraints to matrix
560
  tmpAmat <- rbind(tmpAmat, rbind(c(rep(1, m), rep(0, n+2)),
561
                                  c(rep(1, m), rep(0, n+2))))
562
  
563
  # Add target return to matrix
564
  tmpAmat <- rbind(tmpAmat, as.numeric(targetcon))
565
  
566
  # This step just adds m rows to the matrix to accept box constraints in the next step
567
  tmpAmat <- cbind(tmpAmat, matrix(0, ncol=m, nrow=dim(tmpAmat)[1]))
568
  
569
  # Add lower bound box constraints
570
  tmpAmat <- rbind(tmpAmat, cbind(-diag(m), matrix(0, ncol=n+2, nrow=m), diag(LB)))
571
  
572
  # Add upper bound box constraints
573
  tmpAmat <- rbind(tmpAmat, cbind(diag(m), matrix(0, ncol=n+2, nrow=m), diag(-UB)))
574
  
575
  # Add rows cardinality constraints
576
  tmpAmat <- rbind(tmpAmat, cbind(matrix(0, ncol=m + n + 2, nrow=1), matrix(-1, ncol=m, nrow=1))) 
577
  tmpAmat <- rbind(tmpAmat, cbind(matrix(0, ncol=m + n + 2, nrow=1), matrix(1, ncol=m, nrow=1)))
578
  
579
  # Set up the rhs vector
580
  rhs <- c( rep(0, n), min_sum, max_sum, targetrhs, rep(0, 2*m), -min_pos, max_pos)
581
  
582
  # Set up the dir vector
583
  dir <- c( rep("<=", n), ">=", "<=", targetdir, rep("<=", 2*m), "<=", "<=")
584
  
585
  if(try(!is.null(constraints$groups), silent=TRUE)){
586
    n.groups <- length(constraints$groups)
587
    Amat.group <- matrix(0, nrow=n.groups, ncol=m)
588
    for(i in 1:n.groups){
589
      Amat.group[i, constraints$groups[[i]]] <- 1
590
    }
591
    if(is.null(constraints$cLO)) cLO <- rep(-Inf, n.groups)
592
    if(is.null(constraints$cUP)) cUP <- rep(Inf, n.groups)
593
    zeros <- matrix(0, nrow=n.groups, ncol=(m + n + 2))
594
    tmpAmat <- rbind(tmpAmat, cbind(Amat.group, zeros), cbind(-Amat.group, zeros))
595
    dir <- c(dir, rep(">=", (n.groups + n.groups)))
596
    rhs <- c(rhs, constraints$cLO, -constraints$cUP)
597
  }
598
  
599
  # Add the factor exposures to Amat, dir, and rhs
600
  if(!is.null(constraints$B)){
601
    t.B <- t(constraints$B)
602
    zeros <- matrix(data=0, nrow=nrow(t.B), ncol=(m + n + 2))
603
    tmpAmat <- rbind(tmpAmat, cbind(t.B, zeros), cbind(-t.B, zeros))
604
    dir <- c(dir, rep(">=", 2 * nrow(t.B)))
605
    rhs <- c(rhs, constraints$lower, -constraints$upper)
606
  }
607
  
608
  # Linear objective vector
609
  ROI_objective <- ROI::L_objective(c( rep(0, m), 1, rep(1/n, n) / alpha, 0, rep(0, m)))
610
  
611
  # Set up the types vector with continuous and binary variables
612
  types <- c( rep("C", m), "C", rep("C", n), "C", rep("B", m))
613
  
614
  bnds <- ROI::V_bound( li = 1L:(m + n + 2 + m), lb = c(LB, -1, rep(0, n), 1, rep(0, m)),
615
                   ui = 1L:(m + n + 2 + m), ub = c(UB, 1, rep(Inf, n), 1, rep(1, m)))
616
  
617
  # Set up the optimization problem and solve  
618
  opt.prob <- ROI::OP(objective=ROI_objective, 
619
                 constraints=ROI::L_constraint(L=tmpAmat, dir=dir, rhs=rhs),
620
                 bounds=bnds, types=types)
621
  roi.result <- try(ROI::ROI_solve(x=opt.prob, solver=solver, control=control), silent=TRUE)
622
  if(inherits(roi.result, "try-error")) stop(paste("No solution found:", roi.result))
623
  
624
  # The Rglpk solvers status returns an an integer with status information
625
  # about the solution returned: 0 if the optimal solution was found, a 
626
  #non-zero value otherwise.
627
  if(roi.result$status$code != 0) {
628
    message("Undefined Solution")
629
    return(NULL)
630
  }
631
  
632
  weights <- roi.result$solution[1:m]
633
  names(weights) <- colnames(R)
634
  out <- list()
635
  out$weights <- weights
636
  out$out <- roi.result$objval
637
  es_names <- c("ES", "ETL", "CVaR")
638
  es_idx <- which(es_names %in% names(moments))
639
  obj_vals <- list()
640
  # Calculate the objective values here so that we can use the moments$mean
641
  # and moments$var that might be passed in by the user.
642
  if(!all(moments$mean == 0)){
643
    port.mean <- as.numeric(sum(weights * moments$mean))
644
    names(port.mean) <- "mean"
645
    obj_vals[["mean"]] <- port.mean
646
    port.es <- roi.result$objval
647
    names(port.es) <- es_names[es_idx]
648
    obj_vals[[es_names[es_idx]]] <- port.es
649
  } else {
650
    port.es <- roi.result$objval
651
    names(port.es) <- es_names[es_idx]
652
    obj_vals[[es_names[es_idx]]] <- port.es
653
  }
654
  out$obj_vals <- obj_vals
655
  return(out)
656
}
657

658
##### minimize variance or maximize quadratic utility with turnover constraints #####
659
#' GMV/QU QP Optimization with Turnover Constraint
660
#' 
661
#' This function is called by optimize.portfolio to solve minimum variance or 
662
#' maximum quadratic utility problems with turnover constraint
663
#' 
664
#' @param R xts object of asset returns
665
#' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
666
#' @param moments object of moments computed based on objective functions
667
#' @param lambda risk_aversion parameter
668
#' @param target target return value
669
#' @param init_weights initial weights to compute turnover
670
#' @param solver solver to use
671
#' @param control list of solver control parameters
672
#' @author Ross Bennett
673
gmv_opt_toc <- function(R, constraints, moments, lambda, target, init_weights, solver="quadprog", control=NULL){
674
  # function for minimum variance or max quadratic utility problems
675
  stopifnot("package:corpcor" %in% search() || requireNamespace("corpcor",quietly = TRUE))
676
  stopifnot("package:ROI" %in% search() || requireNamespace("ROI", quietly = TRUE))
677
  plugin <- paste0("ROI.plugin.", solver)
678
  stopifnot(paste0("package:", plugin) %in% search() || requireNamespace(plugin, quietly=TRUE))
679
  
680
  # Check for cleaned returns in moments
681
  if(!is.null(moments$cleanR)) R <- moments$cleanR
682
  
683
  # Modify the returns matrix. This is done because there are 3 sets of
684
  # variables 1) w.initial, 2) w.buy, and 3) w.sell
685
  R0 <- matrix(0, ncol=ncol(R), nrow=nrow(R))
686
  returns <- cbind(R, R0, R0)
687
  V <- cov(returns)
688
  
689
  # number of assets
690
  N <- ncol(R)
691
  
692
  # initial weights for solver
693
  if(is.null(init_weights)) init_weights <- rep(1/ N, N)
694
  
695
  # check for a target return constraint
696
  if(!is.na(target)) {
697
    # If var is the only objective specified, then moments$mean won't be calculated
698
    if(all(moments$mean==0)){
699
      tmp_means <- colMeans(R)
700
    } else {
701
      tmp_means <- moments$mean
702
    }
703
  } else {
704
    tmp_means <- rep(0, N)
705
    target <- 0
706
  }
707
  Amat <- c(tmp_means, rep(0, 2*N))
708
  dir <- "=="
709
  rhs <- target
710
  meq <- N + 1
711
  
712
  # Amat for initial weights
713
  # Amat <- cbind(diag(N), matrix(0, nrow=N, ncol=N*2))
714
  Amat <- rbind(Amat, cbind(diag(N), -1*diag(N), diag(N)))
715
  rhs <- c(rhs, init_weights)
716
  dir <- c(dir, rep("==", N))
717
  
718
  # Amat for turnover constraints
719
  Amat <- rbind(Amat, c(rep(0, N), rep(-1, N), rep(-1, N)))
720
  rhs <- c(rhs, -constraints$turnover_target)
721
  dir <- c(dir, ">=")
722
  
723
  # Amat for positive weights
724
  Amat <- rbind(Amat, cbind(matrix(0, nrow=N, ncol=N), diag(N), matrix(0, nrow=N, ncol=N)))
725
  rhs <- c(rhs, rep(0, N))
726
  dir <- c(dir, rep(">=", N))
727
  
728
  # Amat for negative weights
729
  Amat <- rbind(Amat, cbind(matrix(0, nrow=N, ncol=2*N), diag(N)))
730
  rhs <- c(rhs, rep(0, N))
731
  dir <- c(dir, rep(">=", N))
732
  
733
  # Amat for full investment constraint
734
  Amat <- rbind(Amat, rbind(c(rep(1, N), rep(0,2*N)), c(rep(-1, N), rep(0,2*N))))
735
  rhs <- c(rhs, constraints$min_sum, -constraints$max_sum)
736
  dir <- c(dir, ">=", ">=")
737
  
738
  # Amat for lower box constraints
739
  Amat <- rbind(Amat, cbind(diag(N), diag(0, N), diag(0, N)))
740
  rhs <- c(rhs, constraints$min)
741
  dir <- c(dir, rep(">=", N))
742
  
743
  # Amat for upper box constraints
744
  Amat <- rbind(Amat, cbind(-diag(N), diag(0, N), diag(0, N)))
745
  rhs <- c(rhs, -constraints$max)
746
  dir <- c(dir, rep(">=", N))
747
  
748
  # include group constraints
749
  if(try(!is.null(constraints$groups), silent=TRUE)){
750
    n.groups <- length(constraints$groups)
751
    Amat.group <- matrix(0, nrow=n.groups, ncol=N)
752
    zeros <- matrix(0, nrow=n.groups, ncol=N)
753
    for(i in 1:n.groups){
754
      Amat.group[i, constraints$groups[[i]]] <- 1
755
    }
756
    if(is.null(constraints$cLO)) cLO <- rep(-Inf, n.groups)
757
    if(is.null(constraints$cUP)) cUP <- rep(Inf, n.groups)
758
    Amat <- rbind(Amat, cbind(Amat.group, zeros, zeros))
759
    Amat <- rbind(Amat, cbind(-Amat.group, zeros, zeros))
760
    dir <- c(dir, rep(">=", (n.groups + n.groups)))
761
    rhs <- c(rhs, constraints$cLO, -constraints$cUP)
762
  }
763
  
764
  # Add the factor exposures to Amat, dir, and rhs
765
  if(!is.null(constraints$B)){
766
    t.B <- t(constraints$B)
767
    zeros <- matrix(0, nrow=nrow(t.B), ncol=ncol(t.B))
768
    Amat <- rbind(Amat, cbind(t.B, zeros, zeros))
769
    Amat <- rbind(Amat, cbind(-t.B, zeros, zeros))
770
    dir <- c(dir, rep(">=", 2 * nrow(t.B)))
771
    rhs <- c(rhs, constraints$lower, -constraints$upper)
772
  }
773
  
774
  # Remove the rows of Amat and elements of rhs.vec where rhs is Inf or -Inf
775
  Amat <- Amat[!is.infinite(rhs), ]
776
  rhs <- rhs[!is.infinite(rhs)]
777
  dir <- dir[!is.infinite(rhs)]
778
  
779
  ROI_objective <- ROI::Q_objective(Q=corpcor::make.positive.definite(2*lambda*V), 
780
                               L=rep(-tmp_means, 3))
781
  
782
  opt.prob <- ROI::OP(objective=ROI_objective, 
783
                 constraints=ROI::L_constraint(L=Amat, dir=dir, rhs=rhs))
784
  
785
  roi.result <- try(ROI::ROI_solve(x=opt.prob, solver=solver, control=control), silent=TRUE)
786
  if(inherits(roi.result, "try-error")) stop(paste("No solution found:", roi.result))
787
  
788
  wts <- roi.result$solution
789
  wts.final <- wts[1:N]
790
  
791
  weights <- wts.final
792
  names(weights) <- colnames(R)
793
  out <- list()
794
  out$weights <- weights
795
  out$out <- roi.result$objval
796
  obj_vals <- list()
797
  # Calculate the objective values here so that we can use the moments$mean
798
  # and moments$var that might be passed in by the user.
799
  if(!all(moments$mean == 0)){
800
    port.mean <- as.numeric(sum(weights * moments$mean))
801
    names(port.mean) <- "mean"
802
    obj_vals[["mean"]] <- port.mean
803
    # faster and more efficient way to compute t(w) %*% Sigma %*% w
804
    port.sd <- sqrt(sum(crossprod(weights, moments$var) * weights))
805
    names(port.sd) <- "StdDev"
806
    obj_vals[["StdDev"]] <- port.sd
807
  } else {
808
    # faster and more efficient way to compute t(w) %*% Sigma %*% w
809
    port.sd <- sqrt(sum(crossprod(weights, moments$var) * weights))
810
    names(port.sd) <- "StdDev"
811
    obj_vals[["StdDev"]] <- port.sd
812
  }
813
  out$obj_vals <- obj_vals
814
  return(out)
815
}
816

817
##### minimize variance or maximize quadratic utility with proportional transactioncosts constraints #####
818
#' GMV/QU QP Optimization with Proportional Transaction Cost Constraint
819
#' 
820
#' This function is called by optimize.portfolio to solve minimum variance or 
821
#' maximum quadratic utility problems with proportional transaction cost constraint
822
#' 
823
#' @param R xts object of asset returns
824
#' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
825
#' @param moments object of moments computed based on objective functions
826
#' @param lambda risk_aversion parameter
827
#' @param target target return value
828
#' @param init_weights initial weights to compute turnover
829
#' @param solver solver to use
830
#' @param control list of solver control parameters
831
#' @author Ross Bennett
832
gmv_opt_ptc <- function(R, constraints, moments, lambda, target, init_weights, solver="quadprog", control=NULL){
833
  # function for minimum variance or max quadratic utility problems
834
  # modifying ProportionalCostOpt function from MPO package
835
  stopifnot("package:corpcor" %in% search() || requireNamespace("corpcor", quietly = TRUE))
836
  stopifnot("package:ROI" %in% search() || requireNamespace("ROI", quietly = TRUE))
837
  plugin <- paste0("ROI.plugin.", solver)
838
  stopifnot(paste0("package:", plugin) %in% search() || requireNamespace(plugin, quietly=TRUE))
839
  
840
  # Check for cleaned returns in moments
841
  if(!is.null(moments$cleanR)) R <- moments$cleanR
842
  
843
  # proportional transaction costs
844
  ptc <- constraints$ptc
845
  
846
  # Modify the returns matrix. This is done because there are 3 sets of
847
  # variables 1) w.initial, 2) w.buy, and 3) w.sell
848
  R0 <- matrix(0, ncol=ncol(R), nrow=nrow(R))
849
  returns <- cbind(R, R0, R0)
850
  V <- cov(returns)
851
  
852
  # number of assets
853
  N <- ncol(R)
854
  
855
  # initial weights for solver
856
  if(is.null(init_weights)) init_weights <- rep(1/ N, N)
857
  
858
  # Check for a target return constraint
859
  if(!is.na(target)) {
860
    # If var is the only objective specified, then moments$mean won't be calculated
861
    if(all(moments$mean==0)){
862
      tmp_means <- colMeans(R)
863
    } else {
864
      tmp_means <- moments$mean
865
    }
866
  } else {
867
    tmp_means <- rep(0, N)
868
    target <- 0
869
  }
870
  Amat <- c(tmp_means, rep(0, 2 * N))
871
  dir <- "=="
872
  rhs <- 1 + target
873
  meq <- 1
874
  
875
  # separate the weights into w, w^+, and w^-
876
  # w - w^+ + w^- = 0
877
  Amat <- rbind(Amat, cbind(diag(N), -diag(N), diag(N)))
878
  rhs <- c(rhs, init_weights)
879
  dir <- c(dir, rep("==", N))
880
  meq <- N + 1
881
  
882
  # w+ >= 0
883
  Amat <- rbind(Amat, cbind(diag(0, N), diag(N), diag(0, N)))
884
  rhs <- c(rhs, rep(0, N))
885
  dir <- c(dir, rep(">=", N))
886
  
887
  # w- >= 0
888
  Amat <- rbind(Amat, cbind(diag(0, N), diag(0, N), diag(N)))
889
  rhs <- c(rhs, rep(0, N))
890
  dir <- c(dir, rep(">=", N))
891
  
892
  # 1^T w + tcb^T w^+ + tcs^T w^- >= min_sum
893
  Amat <- rbind(Amat, c(rep(1, N), ptc, ptc))
894
  rhs <- c(rhs, constraints$min_sum)
895
  dir <- c(dir, ">=")
896
  
897
  # 1^T w + tcb^T w^+ + tcs^T w^- <= max_sum
898
  Amat <- rbind(Amat, c(rep(-1, N), -ptc, -ptc))
899
  rhs <- c(rhs, -constraints$max_sum)
900
  dir <- c(dir, ">=")
901
  
902
  # -(1 + tcb)^T w^+ + (1 - tcs)^T w^- >= 0
903
  Amat <- rbind(Amat, c(rep(0, N), -(1 + ptc), (1 - ptc)))
904
  rhs <- c(rhs, 0)
905
  dir <- c(dir, ">=")
906
  
907
  # Amat for lower box constraints
908
  Amat <- rbind(Amat, cbind(diag(N), diag(N), diag(N)))
909
  rhs <- c(rhs, constraints$min)
910
  dir <- c(dir, rep(">=", N))
911
  
912
  # Amat for upper box constraints
913
  Amat <- rbind(Amat, cbind(-diag(N), -diag(N), -diag(N)))
914
  rhs <- c(rhs, -constraints$max)
915
  dir <- c(dir, rep(">=", N))
916
  
917
  # include group constraints
918
  if(try(!is.null(constraints$groups), silent=TRUE)){
919
    n.groups <- length(constraints$groups)
920
    Amat.group <- matrix(0, nrow=n.groups, ncol=N)
921
    for(i in 1:n.groups){
922
      Amat.group[i, constraints$groups[[i]]] <- 1
923
    }
924
    if(is.null(constraints$cLO)) cLO <- rep(-Inf, n.groups)
925
    if(is.null(constraints$cUP)) cUP <- rep(Inf, n.groups)
926
    Amat <- rbind(Amat, cbind(Amat.group, Amat.group, Amat.group))
927
    Amat <- rbind(Amat, cbind(-Amat.group, -Amat.group, -Amat.group))
928
    dir <- c(dir, rep(">=", (n.groups + n.groups)))
929
    rhs <- c(rhs, constraints$cLO, -constraints$cUP)
930
  }
931
  
932
  # Add the factor exposures to Amat, dir, and rhs
933
  if(!is.null(constraints$B)){
934
    t.B <- t(constraints$B)
935
    Amat <- rbind(Amat, cbind(t.B, t.B, t.B))
936
    Amat <- rbind(Amat, cbind(-t.B, -t.B, -t.B))
937
    dir <- c(dir, rep(">=", 2 * nrow(t.B)))
938
    rhs <- c(rhs, constraints$lower, -constraints$upper)
939
  }
940
  d <- c(-tmp_means, rep(0, 2 * N))
941
  
942
  # Remove the rows of Amat and elements of rhs where rhs is Inf or -Inf
943
  Amat <- Amat[!is.infinite(rhs), ]
944
  rhs <- rhs[!is.infinite(rhs)]
945
  dir <- dir[!is.infinite(rhs)]
946
  
947
  ROI_objective <- ROI::Q_objective(Q=corpcor::make.positive.definite(2*lambda*V), L=d)
948
  
949
  opt.prob <- ROI::OP(objective=ROI_objective, 
950
                 constraints=ROI::L_constraint(L=Amat, dir=dir, rhs=rhs))
951
  roi.result <- try(ROI::ROI_solve(x=opt.prob, solver=solver, control=control), silent=TRUE)
952
  if(inherits(roi.result, "try-error")) stop(paste("No solution found:", roi.result))
953
  
954
  wts <- roi.result$solution
955
  weights <- wts[1:N]
956
  names(weights) <- colnames(R)
957
  out <- list()
958
  out$weights <- weights
959
  out$out <- roi.result$objval
960
  obj_vals <- list()
961
  # Calculate the objective values here so that we can use the moments$mean
962
  # and moments$var that might be passed in by the user.
963
  if(!all(moments$mean == 0)){
964
    port.mean <- as.numeric(sum(weights * moments$mean))
965
    names(port.mean) <- "mean"
966
    obj_vals[["mean"]] <- port.mean
967
    # faster and more efficient way to compute t(w) %*% Sigma %*% w
968
    port.sd <- sqrt(sum(crossprod(weights, moments$var) * weights))
969
    names(port.sd) <- "StdDev"
970
    obj_vals[["StdDev"]] <- port.sd
971
  } else {
972
    # faster and more efficient way to compute t(w) %*% Sigma %*% w
973
    port.sd <- sqrt(sum(crossprod(weights, moments$var) * weights))
974
    names(port.sd) <- "StdDev"
975
    obj_vals[["StdDev"]] <- port.sd
976
  }
977
  out$obj_vals <- obj_vals
978
  return(out)
979
}
980

981
##### minimize variance or maximize quadratic utility with leverage constraints #####
982
#' GMV/QU QP Optimization with Turnover Constraint
983
#' 
984
#' This function is called by optimize.portfolio to solve minimum variance or 
985
#' maximum quadratic utility problems with a leverage constraint
986
#' 
987
#' @param R xts object of asset returns
988
#' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
989
#' @param moments object of moments computed based on objective functions
990
#' @param lambda risk_aversion parameter
991
#' @param target target return value
992
#' @param solver solver to use
993
#' @param control list of solver control parameters
994
#' @author Ross Bennett
995
gmv_opt_leverage <- function(R, constraints, moments, lambda, target, solver="quadprog", control=NULL){
996
  # function for minimum variance or max quadratic utility problems
997
  stopifnot("package:corpcor" %in% search() || requireNamespace("corpcor",quietly = TRUE))
998
  stopifnot("package:ROI" %in% search() || requireNamespace("ROI", quietly = TRUE))
999
  plugin <- paste0("ROI.plugin.", solver)
1000
  stopifnot(paste0("package:", plugin) %in% search() || requireNamespace(plugin, quietly=TRUE))
1001
  
1002
  # Check for cleaned returns in moments
1003
  if(!is.null(moments$cleanR)) R <- moments$cleanR
1004
  
1005
  # Modify the returns matrix. This is done because there are 3 sets of
1006
  # variables 1) w.initial, 2) w.buy, and 3) w.sell
1007
  R0 <- matrix(0, ncol=ncol(R), nrow=nrow(R))
1008
  returns <- cbind(R, R0, R0)
1009
  V <- cov(returns)
1010
  
1011
  # number of assets
1012
  N <- ncol(R)
1013
  
1014
  # check for a target return constraint
1015
  if(!is.na(target)) {
1016
    # If var is the only objective specified, then moments$mean won't be calculated
1017
    if(all(moments$mean==0)){
1018
      tmp_means <- colMeans(R)
1019
    } else {
1020
      tmp_means <- moments$mean
1021
    }
1022
  } else {
1023
    tmp_means <- rep(0, N)
1024
    target <- 0
1025
  }
1026
  Amat <- c(tmp_means, rep(0, 2*N))
1027
  dir <- "=="
1028
  rhs <- target
1029
  # meq <- N + 1
1030
  
1031
  # separate the weights into w, w+, and w-
1032
  # w - w+ + w- = 0
1033
  Amat <- rbind(Amat, cbind(diag(N), -1*diag(N), diag(N)))
1034
  rhs <- c(rhs, rep(0, N))
1035
  dir <- c(dir, rep("==", N))
1036
  
1037
  # Amat for leverage constraints
1038
  Amat <- rbind(Amat, c(rep(0, N), rep(-1, N), rep(-1, N)))
1039
  rhs <- c(rhs, -constraints$leverage)
1040
  dir <- c(dir, ">=")
1041
  
1042
  # Amat for positive weights
1043
  Amat <- rbind(Amat, cbind(matrix(0, nrow=N, ncol=N), diag(N), matrix(0, nrow=N, ncol=N)))
1044
  rhs <- c(rhs, rep(0, N))
1045
  dir <- c(dir, rep(">=", N))
1046
  
1047
  # Amat for negative weights
1048
  Amat <- rbind(Amat, cbind(matrix(0, nrow=N, ncol=2*N), diag(N)))
1049
  rhs <- c(rhs, rep(0, N))
1050
  dir <- c(dir, rep(">=", N))
1051
  
1052
  # Amat for full investment constraint
1053
  Amat <- rbind(Amat, rbind(c(rep(1, N), rep(0,2*N)), 
1054
                            c(rep(-1, N), rep(0,2*N))))
1055
  rhs <- c(rhs, constraints$min_sum, -constraints$max_sum)
1056
  dir <- c(dir, ">=", ">=")
1057
  
1058
  # Amat for lower box constraints
1059
  Amat <- rbind(Amat, cbind(diag(N), diag(0, N), diag(0, N)))
1060
  rhs <- c(rhs, constraints$min)
1061
  dir <- c(dir, rep(">=", N))
1062
  
1063
  # Amat for upper box constraints
1064
  Amat <- rbind(Amat, cbind(-diag(N), diag(0, N), diag(0, N)))
1065
  rhs <- c(rhs, -constraints$max)
1066
  dir <- c(dir, rep(">=", N))
1067
  
1068
  # include group constraints
1069
  if(try(!is.null(constraints$groups), silent=TRUE)){
1070
    n.groups <- length(constraints$groups)
1071
    Amat.group <- matrix(0, nrow=n.groups, ncol=N)
1072
    zeros <- matrix(0, nrow=n.groups, ncol=N)
1073
    for(i in 1:n.groups){
1074
      Amat.group[i, constraints$groups[[i]]] <- 1
1075
    }
1076
    if(is.null(constraints$cLO)) cLO <- rep(-Inf, n.groups)
1077
    if(is.null(constraints$cUP)) cUP <- rep(Inf, n.groups)
1078
    Amat <- rbind(Amat, cbind(Amat.group, zeros, zeros))
1079
    Amat <- rbind(Amat, cbind(-Amat.group, zeros, zeros))
1080
    dir <- c(dir, rep(">=", (n.groups + n.groups)))
1081
    rhs <- c(rhs, constraints$cLO, -constraints$cUP)
1082
  }
1083
  
1084
  # Add the factor exposures to Amat, dir, and rhs
1085
  if(!is.null(constraints$B)){
1086
    t.B <- t(constraints$B)
1087
    zeros <- matrix(0, nrow=nrow(t.B), ncol=ncol(t.B))
1088
    Amat <- rbind(Amat, cbind(t.B, zeros, zeros))
1089
    Amat <- rbind(Amat, cbind(-t.B, zeros, zeros))
1090
    dir <- c(dir, rep(">=", 2 * nrow(t.B)))
1091
    rhs <- c(rhs, constraints$lower, -constraints$upper)
1092
  }
1093
  
1094
  # Remove the rows of Amat and elements of rhs.vec where rhs is Inf or -Inf
1095
  Amat <- Amat[!is.infinite(rhs), ]
1096
  rhs <- rhs[!is.infinite(rhs)]
1097
  dir <- dir[!is.infinite(rhs)]
1098
  
1099
  ROI_objective <- ROI::Q_objective(Q=corpcor::make.positive.definite(2*lambda*V), 
1100
                               L=rep(-tmp_means, 3))
1101
  
1102
  opt.prob <- ROI::OP(objective=ROI_objective, 
1103
                 constraints=ROI::L_constraint(L=Amat, dir=dir, rhs=rhs))
1104
  
1105
  roi.result <- try(ROI::ROI_solve(x=opt.prob, solver=solver, control=control), silent=TRUE)
1106
  if(inherits(roi.result, "try-error")) stop(paste("No solution found:", roi.result))
1107
  
1108
  wts <- roi.result$solution
1109
  wts.final <- wts[1:N]
1110
  
1111
  weights <- wts.final
1112
  names(weights) <- colnames(R)
1113
  out <- list()
1114
  out$weights <- weights
1115
  out$out <- roi.result$objval
1116
  obj_vals <- list()
1117
  # Calculate the objective values here so that we can use the moments$mean
1118
  # and moments$var that might be passed in by the user.
1119
  if(!all(moments$mean == 0)){
1120
    port.mean <- as.numeric(sum(weights * moments$mean))
1121
    names(port.mean) <- "mean"
1122
    obj_vals[["mean"]] <- port.mean
1123
    # faster and more efficient way to compute t(w) %*% Sigma %*% w
1124
    port.sd <- sqrt(sum(crossprod(weights, moments$var) * weights))
1125
    names(port.sd) <- "StdDev"
1126
    obj_vals[["StdDev"]] <- port.sd
1127
  } else {
1128
    # faster and more efficient way to compute t(w) %*% Sigma %*% w
1129
    port.sd <- sqrt(sum(crossprod(weights, moments$var) * weights))
1130
    names(port.sd) <- "StdDev"
1131
    obj_vals[["StdDev"]] <- port.sd
1132
  }
1133
  out$obj_vals <- obj_vals
1134
  return(out)
1135
}
1136

1137
# This function uses optimize() to find the target return value that 
1138
# results in the maximum starr ratio (mean / ES).
1139
# returns the target return value
1140
mean_etl_opt <- function(R, constraints, moments, alpha, solver, control){
1141
  # create a copy of the moments that can be modified
1142
  tmp_moments <- moments
1143
  
1144
  # Find the maximum return
1145
  if(!is.null(constraints$max_pos)){
1146
    max_ret <- maxret_milp_opt(R=R, constraints=constraints, moments=moments, target=NA, solver=solver, control=control)
1147
  } else {
1148
    max_ret <- maxret_opt(R=R, moments=moments, constraints=constraints, target=NA, solver=solver, control=control)
1149
  }
1150
  max_mean <- as.numeric(-max_ret$out)
1151
  
1152
  # Find the minimum return
1153
  tmp_moments$mean <- -1 * moments$mean
1154
  if(!is.null(constraints$max_pos)){
1155
    min_ret <- maxret_milp_opt(R=R, constraints=constraints, moments=tmp_moments, target=NA, solver=solver, control=control)
1156
  } else {
1157
    min_ret <- maxret_opt(R=R, constraints=constraints, moments=tmp_moments, target=NA, solver=solver, control=control)
1158
  }
1159
  min_mean <- as.numeric(min_ret$out)
1160
  
1161
  # use optimize() to find the target return value that maximizes sharpe ratio
1162
  opt <- try(optimize(f=starr_obj_fun, R=R, constraints=constraints, 
1163
                      solver=solver, moments=moments, alpha=alpha,
1164
                      lower=min_mean, upper=max_mean, control=control,
1165
                      maximum=TRUE, tol=.Machine$double.eps), 
1166
             silent=TRUE)
1167
  if(inherits(opt, "try-error")){
1168
    stop(paste("Objective function failed with message\n", opt))
1169
    return(NULL)
1170
  }
1171
  return(opt$maximum)
1172
}
1173

1174
# Function to calculate the starr ratio.
1175
# Used as the objective function for optimize()
1176
starr_obj_fun <- function(target_return, R, constraints, moments, alpha, solver, control){
1177
  if(!is.null(constraints$max_pos)){
1178
    opt <- etl_milp_opt(R=R, constraints=constraints, moments=moments, 
1179
                        target=target_return, alpha=alpha, solver=solver, 
1180
                        control=control)
1181
  } else {
1182
    opt <- etl_opt(R=R, constraints=constraints, moments=moments, 
1183
                   target=target_return, alpha=alpha, solver=solver,
1184
                   control=control)
1185
  }
1186
  weights <- matrix(opt$weights, ncol=1)
1187
  opt_mean <- sum(weights * moments$mean)
1188
  opt_etl <- as.numeric(opt$out)
1189
  starr <- opt_mean / opt_etl
1190
  return(starr)
1191
}
1192

1193

1194
# This was my implementation of a binary search for the maximum starr ratio 
1195
# target return. Better to use optimize() in R rather than my method. -Ross Bennett
1196
# mean_etl_opt <- function(R, constraints, moments, target, alpha, solver="glpk", tol=.Machine$double.eps^0.5, maxit=50){
1197
#   # This function returns the target mean return that maximizes mean / etl (i.e. starr)
1198
#   
1199
#   # if all(moments$mean == 0) then the user did not specify mean as an objective,
1200
#   # and we just want to return the target mean return value
1201
#   if(all(moments$mean == 0)) return(target)
1202
#   
1203
#   fmean <- matrix(moments$mean, ncol=1)
1204
#   
1205
#   # can't use optimize.portfolio here, this function is called inside 
1206
#   # optimize.portfolio and will throw an error message about nesting too deeply
1207
#   
1208
#   # Find the maximum return
1209
#   if(!is.null(constraints$max_pos)){
1210
#     max_ret <- maxret_milp_opt(R=R, constraints=constraints, moments=moments, target=NA, solver=solver)
1211
#   } else {
1212
#     max_ret <- maxret_opt(R=R, moments=moments, constraints=constraints, target=NA, solver=solver)
1213
#   }
1214
#   max_mean <- as.numeric(-max_ret$out)
1215
#   
1216
#   # Find the starr at the maximum etl portfolio
1217
#   if(!is.null(constraints$max_pos)){
1218
#     ub_etl <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=max_mean, alpha=alpha, solver=solver)
1219
#   } else {
1220
#     ub_etl <- etl_opt(R=R, constraints=constraints, moments=moments, target=max_mean, alpha=alpha, solver=solver)
1221
#   }
1222
#   ub_weights <- matrix(ub_etl$weights, ncol=1)
1223
#   ub_mean <- as.numeric(t(ub_weights) %*% fmean)
1224
#   ub_etl <- as.numeric(ub_etl$out)
1225
#   # starr at the upper bound
1226
#   ub_starr <- ub_mean / ub_etl
1227
#   if(is.infinite(ub_starr)) stop("Inf value for STARR, objective value is 0")
1228
#   
1229
#   # cat("ub_mean", ub_mean, "\n")
1230
#   # cat("ub_etl", ub_etl, "\n")
1231
#   # cat("ub_starr", ub_starr, "\n")
1232
#   
1233
#   # Find the starr at the minimum etl portfolio
1234
#   if(!is.null(constraints$max_pos)){
1235
#     lb_etl <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=NA, alpha=alpha, solver=solver)
1236
#   } else {
1237
#     lb_etl <- etl_opt(R=R, constraints=constraints, moments=moments, target=NA, alpha=alpha, solver=solver)
1238
#   }
1239
#   lb_weights <- matrix(lb_etl$weights)  
1240
#   lb_mean <- as.numeric(t(lb_weights) %*% fmean)  
1241
#   lb_etl <- as.numeric(lb_etl$out)
1242
#   
1243
#   # starr at the lower bound
1244
#   lb_starr <- lb_mean / lb_etl
1245
#   # if(is.infinite(lb_starr)) stop("Inf value for STARR, objective value is 0")
1246
#   
1247
#   # set lb_starr equal to 0, should this be a negative number like -1e6?
1248
#   # the lb_* values will be 0 for a dollar-neutral strategy so we need to reset the values
1249
#   if(is.na(lb_starr) | is.infinite(lb_starr)) lb_starr <- 0
1250
#   
1251
#   # cat("lb_mean", lb_mean, "\n")
1252
#   # cat("lb_etl", lb_etl, "\n")
1253
#   # cat("lb_starr", lb_starr, "\n")
1254
#   
1255
#   # want to find the return that maximizes mean / etl
1256
#   i <- 1
1257
#   while((abs(ub_starr - lb_starr) > tol) & (i < maxit)){
1258
#     # bisection method to find the maximum mean / etl
1259
#     
1260
#     # print(i)
1261
#     # cat("ub_starr", ub_starr, "\n")
1262
#     # cat("lb_starr", lb_starr, "\n")
1263
#     # print("**********")
1264
#     # Find the starr at the mean return midpoint
1265
#     new_ret <- (lb_mean + ub_mean) / 2
1266
#     if(!is.null(constraints$max_pos)){
1267
#       mid <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha, solver=solver)
1268
#     } else {
1269
#       mid <- etl_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha, solver=solver)
1270
#     }
1271
#     # print(mid)
1272
#     mid_weights <- matrix(mid$weights, ncol=1)
1273
#     mid_mean <- as.numeric(t(mid_weights) %*% fmean)
1274
#     mid_etl <- as.numeric(mid$out)
1275
#     mid_starr <- mid_mean / mid_etl
1276
#     # the mid_* values MIGHT be 0 for a dollar-neutral strategy so we need to reset the values
1277
#     # if(is.na(mid_starr) | is.infinite(mid_starr)) mid_starr <- 0
1278
#     # tmp_starr <- mid_starr
1279
#     
1280
#     # cat("mid_mean", mid_mean, "\n")
1281
#     # cat("mid_etl", mid_etl, "\n")
1282
#     # cat("mid_starr", mid_starr, "\n")
1283
#     
1284
#     if(mid_starr > ub_starr){
1285
#       # if mid_starr > ub_starr then mid_starr becomes the new upper bound
1286
#       ub_mean <- mid_mean
1287
#       ub_starr <- mid_starr
1288
#       new_ret <- (lb_mean + ub_mean) / 2
1289
#       if(!is.null(constraints$max_pos)){
1290
#         mid <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha, solver=solver)
1291
#       } else {
1292
#         mid <- etl_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha, solver=solver)
1293
#       }
1294
#       mid_weights <- matrix(mid$weights, ncol=1)
1295
#       mid_mean <- as.numeric(t(mid_weights) %*% fmean)
1296
#       mid_etl <- as.numeric(mid$out)
1297
#       mid_starr <- mid_mean / mid_etl
1298
#       # the mid_* values MIGHT be 0 for a dollar-neutral strategy so we need to reset the values
1299
#       # if(is.na(mid_starr) | is.infinite(mid_starr)) mid_starr <- 0
1300
#     } else if(mid_starr > lb_starr){
1301
#       # if mid_starr > lb_starr then mid_starr becomes the new lower bound
1302
#       lb_mean <- mid_mean
1303
#       lb_starr <- mid_starr
1304
#       new_ret <- (lb_mean + ub_mean) / 2
1305
#       if(!is.null(constraints$max_pos)){
1306
#         mid <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha, solver=solver)
1307
#       } else {
1308
#         mid <- etl_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha, solver=solver)
1309
#       }
1310
#       mid_weights <- matrix(mid$weights, ncol=1)
1311
#       mid_mean <- as.numeric(t(mid_weights) %*% fmean)
1312
#       mid_etl <- as.numeric(mid$out)
1313
#       mid_starr <- mid_mean / mid_etl
1314
#       # the mid_* values MIGHT be 0 for a dollar-neutral strategy so we need to reset the values
1315
#       # if(is.na(mid_starr) | is.infinite(mid_starr)) mid_starr <- 0
1316
#     }
1317
#     i <- i + 1
1318
#   }
1319
#   return(new_ret)
1320
# }
1321

1322

1323

1324
# Function to calculate the sharpe ratio.
1325
# Used as the objective function for optimize()
1326
sharpe_obj_fun <- function(target_return, R, constraints, moments, lambda_hhi=NULL, conc_groups=NULL, solver="quadprog", control=control){
1327
  opt <- gmv_opt(R=R, constraints=constraints, moments=moments, lambda=1, 
1328
                 target=target_return, lambda_hhi=lambda_hhi,
1329
                 conc_groups=conc_groups, solver=solver, control=control)
1330
  weights <- opt$weights
1331
  # opt_mean <- as.numeric(t(weights) %*% matrix(moments$mean, ncol=1))
1332
  opt_mean <- sum(weights * moments$mean)
1333
  # opt_sd <- as.numeric(sqrt(t(weights) %*% moments$var %*% weights))
1334
  # faster and more efficient way to compute t(w) %*% Sigma %*% w
1335
  opt_sd <- sqrt(sum(crossprod(weights, moments$var) * weights))
1336
  opt_sr <- opt_mean / opt_sd
1337
  return(opt_sr)
1338
}
1339

1340
# This function uses optimize() to find the target return value that 
1341
# results in the maximum sharpe ratio (mean / sd).
1342
# returns the target return value
1343
max_sr_opt <- function(R, constraints, moments, lambda_hhi, conc_groups, solver, control){
1344
  # create a copy of the moments that can be modified
1345
  tmp_moments <- moments
1346
  
1347
  # Find the maximum return
1348
  max_ret <- maxret_opt(R=R, moments=moments, constraints=constraints, 
1349
                        target=NA, solver="glpk", control=control)
1350
  max_mean <- as.numeric(-max_ret$out)
1351
  
1352
  # Find the minimum return
1353
  tmp_moments$mean <- -1 * moments$mean
1354
  min_ret <- maxret_opt(R=R, moments=tmp_moments, constraints=constraints, 
1355
                        target=NA, solver="glpk", control=control)
1356
  min_mean <- as.numeric(min_ret$out)
1357
  
1358
  # use optimize() to find the target return value that maximizes sharpe ratio
1359
  opt <- try(optimize(f=sharpe_obj_fun, R=R, constraints=constraints, 
1360
                      solver=solver, lambda_hhi=lambda_hhi, 
1361
                      conc_groups=conc_groups, moments=moments, control=control,
1362
                      lower=min_mean, upper=max_mean, 
1363
                      maximum=TRUE, tol=.Machine$double.eps), 
1364
             silent=TRUE)
1365
  if(inherits(opt, "try-error")){
1366
    stop(paste("Objective function failed with message\n", opt))
1367
    return(NULL)
1368
  }
1369
  return(opt$maximum)
1370
}
1371

1372

1373
# This was my implementation of a binary search for the maximum sharpe ratio 
1374
# target return. Better to use optimize() in R rather than my method. -Ross Bennett
1375
# max_sr_opt <- function(R, constraints, moments, lambda, target, lambda_hhi, conc_groups, solver="quadprog", tol=.Machine$double.eps^0.5, maxit=50){
1376
#   # This function returns the target mean return that maximizes mean / sd (i.e. sharpe ratio)
1377
#   
1378
#   # get the forecast mean from moments
1379
#   fmean <- matrix(moments$mean, ncol=1)
1380
#   
1381
#   # Find the maximum return
1382
#   max_ret <- maxret_opt(R=R, moments=moments, constraints=constraints, target=NA)
1383
#   max_mean <- as.numeric(-max_ret$out)
1384
#   
1385
#   # Calculate the sr at the maximum mean return portfolio
1386
#   ub_weights <- matrix(max_ret$weights, ncol=1)
1387
#   ub_mean <- max_mean
1388
#   ub_sd <- as.numeric(sqrt(t(ub_weights) %*% moments$var %*% ub_weights))
1389
#   # sr at the upper bound
1390
#   ub_sr <- ub_mean / ub_sd
1391
#   
1392
#   # Calculate the sr at the miminum var portfolio
1393
#   tmpmoments <- moments
1394
#   tmpmoments$mean <- rep(0, length(moments$mean))
1395
#   lb_sr <- gmv_opt(R=R, constraints=constraints, moments=tmpmoments, lambda=1, target=NA, lambda_hhi=lambda_hhi, conc_groups=conc_groups, solver=solver)
1396
#   lb_weights <- matrix(lb_sr$weights)
1397
#   lb_mean <- as.numeric(t(lb_weights) %*% fmean)
1398
#   lb_sd <- as.numeric(sqrt(t(lb_weights) %*% moments$var %*% lb_weights))
1399
#   # sr at the lower bound
1400
#   lb_sr <- lb_mean / lb_sd
1401
#   
1402
#   # cat("lb_mean:", lb_mean, "\n")
1403
#   # cat("ub_mean:", ub_mean, "\n")
1404
#   # print("**********")
1405
#   
1406
#   # want to find the return that maximizes mean / sd
1407
#   i <- 1
1408
#   while((abs(ub_sr - lb_sr) > tol) & (i < maxit)){
1409
#     # bisection method to find the maximum mean / sd
1410
#     
1411
#     # Find the starr at the mean return midpoint
1412
#     new_ret <- (lb_mean + ub_mean) / 2
1413
#     mid <- gmv_opt(R=R, constraints=constraints, moments=tmpmoments, lambda=1, target=new_ret, lambda_hhi=lambda_hhi, conc_groups=conc_groups, solver=solver)
1414
#     mid_weights <- matrix(mid$weights, ncol=1)
1415
#     mid_mean <- as.numeric(t(mid_weights) %*% fmean)
1416
#     mid_sd <- as.numeric(sqrt(t(mid_weights) %*% moments$var %*% mid_weights))
1417
#     mid_sr <- mid_mean / mid_sd
1418
#     # tmp_sr <- mid_sr
1419
#     
1420
#     # print(i)
1421
#     # cat("new_ret:", new_ret, "\n")
1422
#     # cat("mid_sr:", mid_sr, "\n")
1423
#     # print("**********")
1424
#     
1425
#     if(mid_sr > ub_sr){
1426
#       # if mid_sr > ub_sr then mid_sr becomes the new upper bound
1427
#       ub_mean <- mid_mean
1428
#       ub_sr <- mid_sr
1429
#       new_ret <- (lb_mean + ub_mean) / 2
1430
#       mid <- gmv_opt(R=R, constraints=constraints, moments=tmpmoments, lambda=1, target=new_ret, lambda_hhi=lambda_hhi, conc_groups=conc_groups, solver=solver)
1431
#       mid_weights <- matrix(mid$weights, ncol=1)
1432
#       mid_mean <- as.numeric(t(mid_weights) %*% fmean)
1433
#       mid_sd <- as.numeric(sqrt(t(mid_weights) %*% moments$var %*% mid_weights))
1434
#       mid_sr <- mid_mean / mid_sd
1435
#     } else if(mid_sr > lb_sr){
1436
#       # if mid_sr > lb_sr then mid_sr becomes the new lower bound
1437
#       lb_mean <- mid_mean
1438
#       lb_sr <- mid_sr
1439
#       new_ret <- (lb_mean + ub_mean) / 2
1440
#       mid <- gmv_opt(R=R, constraints=constraints, moments=tmpmoments, lambda=1, target=new_ret, lambda_hhi=lambda_hhi, conc_groups=conc_groups, solver=solver)
1441
#       mid_weights <- matrix(mid$weights, ncol=1)
1442
#       mid_mean <- as.numeric(t(mid_weights) %*% fmean)
1443
#       mid_sd <- as.numeric(sqrt(t(mid_weights) %*% moments$var %*% mid_weights))
1444
#       mid_sr <- mid_mean / mid_sd
1445
#     }
1446
#     i <- i + 1
1447
#   }
1448
#   return(new_ret)
1449
# }
1450

1451

1452
###############################################################################
1453
# R (https://r-project.org/) Numeric Methods for Optimization of Portfolios
1454
#
1455
# Copyright (c) 2004-2021 Brian G. Peterson, Peter Carl, Ross Bennett, Kris Boudt
1456
#
1457
# This library is distributed under the terms of the GNU Public License (GPL)
1458
# for full details see the file COPYING
1459
#
1460
# $Id$
1461
#
1462
###############################################################################
1463

1464
Product

Resources

Company
