##### GMV and QU QP Function #####
#' GMV/QU QP Optimization
#'
#' This function is called by optimize.portfolio to solve minimum variance or
#' maximum quadratic utility problems
#'
#' @param R xts object of asset returns
#' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
#' @param moments object of moments computed based on objective functions
#' @param lambda risk_aversion parameter
#' @param target target return value
#' @param lambda_hhi concentration aversion parameter
#' @param conc_groups list of vectors specifying the groups of the assets.
#' @param solver solver to use
#' @param control list of solver control parameters
#' @author Ross Bennett
gmv_opt <- function(R, constraints, moments, lambda, target, lambda_hhi, conc_groups, solver="quadprog", control=NULL){
stopifnot("package:ROI" %in% search() || require("ROI", quietly=TRUE))
plugin <- paste0("ROI.plugin.", solver)
stopifnot(paste0("package:", plugin) %in% search() || require(plugin, quietly=TRUE, character.only=TRUE))
# Check for cleaned returns in moments
if(!is.null(moments$cleanR)) R <- moments$cleanR
# Number of assets
N <- ncol(R)
# Check for a target return constraint
if(!is.na(target)) {
# If var is the only objective specified, then moments$mean won't be calculated
if(all(moments$mean==0)){
tmp_means <- colMeans(R)
} else {
tmp_means <- moments$mean
}
} else {
tmp_means <- rep(0, N)
target <- 0
}
Amat <- tmp_means
dir.vec <- "=="
rhs.vec <- target
meq <- 1
# Set up initial A matrix for leverage constraints
Amat <- rbind(Amat, rep(1, N), rep(-1, N))
dir.vec <- c(dir.vec, ">=",">=")
rhs.vec <- c(rhs.vec, constraints$min_sum, -constraints$max_sum)
# Add min box constraints
Amat <- rbind(Amat, diag(N))
dir.vec <- c(dir.vec, rep(">=", N))
rhs.vec <- c(rhs.vec, constraints$min)
# Add max box constraints
Amat <- rbind(Amat, -1*diag(N))
dir.vec <- c(dir.vec, rep(">=", N))
rhs.vec <- c(rhs.vec, -constraints$max)
# Applying box constraints
lb <- constraints$min
ub <- constraints$max
bnds <- V_bound(li=seq.int(1L, N), lb=as.numeric(lb),
ui=seq.int(1L, N), ub=as.numeric(ub))
# Include group constraints
if(try(!is.null(constraints$groups), silent=TRUE)){
n.groups <- length(constraints$groups)
Amat.group <- matrix(0, nrow=n.groups, ncol=N)
for(i in 1:n.groups){
Amat.group[i, constraints$groups[[i]]] <- 1
}
if(is.null(constraints$cLO)) cLO <- rep(-Inf, n.groups)
if(is.null(constraints$cUP)) cUP <- rep(Inf, n.groups)
Amat <- rbind(Amat, Amat.group, -Amat.group)
dir.vec <- c(dir.vec, rep(">=", (n.groups + n.groups)))
rhs.vec <- c(rhs.vec, constraints$cLO, -constraints$cUP)
}
# Add the factor exposures to Amat, dir.vec, and rhs.vec
if(!is.null(constraints$B)){
t.B <- t(constraints$B)
Amat <- rbind(Amat, t.B, -t.B)
dir.vec <- c(dir.vec, rep(">=", 2 * nrow(t.B)))
rhs.vec <- c(rhs.vec, constraints$lower, -constraints$upper)
}
# quadprog cannot handle infinite values so replace Inf with .Machine$double.xmax
# This is the strategy used in ROI
# Amat[ is.infinite(Amat) & (Amat <= 0) ] <- -.Machine$double.xmax
# Amat[ is.infinite(Amat) & (Amat >= 0) ] <- .Machine$double.xmax
# rhs.vec[is.infinite(rhs.vec) & (rhs.vec <= 0)] <- -.Machine$double.xmax
# rhs.vec[is.infinite(rhs.vec) & (rhs.vec >= 0)] <- .Machine$double.xmax
# Remove the rows of Amat and elements of dir.vec and rhs.vec where rhs.vec is Inf or -Inf
Amat <- Amat[!is.infinite(rhs.vec), ]
dir.vec <- dir.vec[!is.infinite(rhs.vec)]
rhs.vec <- rhs.vec[!is.infinite(rhs.vec)]
# Set up the quadratic objective
if(!is.null(lambda_hhi)){
if(length(lambda_hhi) == 1 & is.null(conc_groups)){
ROI_objective <- Q_objective(Q=2*lambda*(moments$var + lambda_hhi * diag(N)), L=-moments$mean) # ROI
#Dmat <- 2*lambda*(moments$var + lambda_hhi * diag(N)) # solve.QP
#dvec <- moments$mean # solve.QP
} else if(!is.null(conc_groups)){
# construct the matrix with concentration aversion values by group
hhi_mat <- matrix(0, nrow=N, ncol=N)
vec <- 1:N
for(i in 1:length(conc_groups)){
tmpI <- diag(N)
tmpvec <- conc_groups[[i]]
zerovec <- setdiff(vec, tmpvec)
for(j in 1:length(zerovec)){
tmpI[zerovec[j], ] <- rep(0, N)
}
hhi_mat <- hhi_mat + lambda_hhi[i] * tmpI
}
ROI_objective <- Q_objective(Q=2*lambda*(moments$var + hhi_mat), L=-moments$mean) # ROI
#Dmat <- 2 * lambda * (moments$var + hhi_mat) # solve.QP
#dvec <- moments$mean # solve.QP
}
} else {
ROI_objective <- Q_objective(Q=2*lambda*moments$var, L=-moments$mean) # ROI
#Dmat <- 2 * lambda * moments$var # solve.QP
#dvec <- moments$mean # solve.QP
}
# set up the optimization problem and solve
opt.prob <- OP(objective=ROI_objective,
constraints=L_constraint(L=Amat, dir=dir.vec, rhs=rhs.vec),
bounds=bnds)
result <- try(ROI_solve(x=opt.prob, solver=solver, control=control), silent=TRUE)
# result <- try(solve.QP(Dmat=Dmat, dvec=dvec, Amat=t(Amat), bvec=rhs.vec, meq=meq), silent=TRUE)
if(inherits(x=result, "try-error")) stop(paste("No solution found:", result))
weights <- result$solution[1:N]
names(weights) <- colnames(R)
out <- list()
out$weights <- weights
out$out <- result$objval
obj_vals <- list()
# Calculate the objective values here so that we can use the moments$mean
# and moments$var that might be passed in by the user. This will avoid
# the extra call to constrained_objective
if(!all(moments$mean == 0)){
port.mean <- as.numeric(sum(weights * moments$mean))
names(port.mean) <- "mean"
obj_vals[["mean"]] <- port.mean
# faster and more efficient way to compute t(w) %*% Sigma %*% w
port.sd <- sqrt(sum(crossprod(weights, moments$var) * weights))
names(port.sd) <- "StdDev"
obj_vals[["StdDev"]] <- port.sd
} else {
# faster and more efficient way to compute t(w) %*% Sigma %*% w
port.sd <- sqrt(sum(crossprod(weights, moments$var) * weights))
names(port.sd) <- "StdDev"
obj_vals[["StdDev"]] <- port.sd
}
out$obj_vals <- obj_vals
# out$out <- result$objval # ROI
# out$call <- call # need to get the call outside of the function
return(out)
}
##### Maximize Return LP Function #####
#' Maximum Return LP Optimization
#'
#' This function is called by optimize.portfolio to solve maximum return
#'
#' @param R xts object of asset returns
#' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
#' @param moments object of moments computed based on objective functions
#' @param target target return value
#' @param solver solver to use
#' @param control list of solver control parameters
#' @author Ross Bennett
maxret_opt <- function(R, moments, constraints, target, solver="glpk", control=NULL){
stopifnot("package:ROI" %in% search() || require("ROI",quietly = TRUE))
plugin <- paste0("ROI.plugin.", solver)
stopifnot(paste0("package:", plugin) %in% search() || require(plugin, quietly=TRUE, character.only=TRUE))
# Check for cleaned returns in moments
if(!is.null(moments$cleanR)) R <- moments$cleanR
N <- ncol(R)
# Applying box constraints
# maxret_opt needs non infinite values for upper and lower bounds
lb <- constraints$min
ub <- constraints$max
if(any(is.infinite(lb)) | any(is.infinite(ub))){
warning("Inf or -Inf values detected in box constraints, maximum return
objectives must have finite box constraint values.")
ub[is.infinite(ub)] <- max(abs(c(constraints$min_sum, constraints$max_sum)))
lb[is.infinite(lb)] <- 0
}
bnds <- V_bound(li=seq.int(1L, N), lb=as.numeric(lb),
ui=seq.int(1L, N), ub=as.numeric(ub))
# set up initial A matrix for leverage constraints
Amat <- rbind(rep(1, N), rep(1, N))
dir.vec <- c(">=","<=")
rhs.vec <- c(constraints$min_sum, constraints$max_sum)
# check for a target return
if(!is.na(target)) {
Amat <- rbind(Amat, moments$mean)
dir.vec <- c(dir.vec, "==")
rhs.vec <- c(rhs.vec, target)
}
# include group constraints
if(try(!is.null(constraints$groups), silent=TRUE)){
n.groups <- length(constraints$groups)
Amat.group <- matrix(0, nrow=n.groups, ncol=N)
for(i in 1:n.groups){
Amat.group[i, constraints$groups[[i]]] <- 1
}
if(is.null(constraints$cLO)) cLO <- rep(-Inf, n.groups)
if(is.null(constraints$cUP)) cUP <- rep(Inf, n.groups)
Amat <- rbind(Amat, Amat.group, -Amat.group)
dir.vec <- c(dir.vec, rep(">=", (n.groups + n.groups)))
rhs.vec <- c(rhs.vec, constraints$cLO, -constraints$cUP)
}
# Add the factor exposures to Amat, dir.vec, and rhs.vec
if(!is.null(constraints$B)){
t.B <- t(constraints$B)
Amat <- rbind(Amat, t.B, -t.B)
dir.vec <- c(dir.vec, rep(">=", 2 * nrow(t.B)))
rhs.vec <- c(rhs.vec, constraints$lower, -constraints$upper)
}
# set up the linear objective
ROI_objective <- L_objective(L=-moments$mean)
# objL <- -moments$mean
# set up the optimization problem and solve
opt.prob <- OP(objective=ROI_objective,
constraints=L_constraint(L=Amat, dir=dir.vec, rhs=rhs.vec),
bounds=bnds)
roi.result <- try(ROI_solve(x=opt.prob, solver=solver, control=control), silent=TRUE)
if(inherits(roi.result, "try-error")) stop(paste("No solution found:", roi.result))
# roi.result <- Rglpk_solve_LP(obj=objL, mat=Amat, dir=dir.vec, rhs=rhs.vec, bounds=bnds)
# The Rglpk solvers status returns an an integer with status information
# about the solution returned: 0 if the optimal solution was found, a
#non-zero value otherwise.
if(roi.result$status$code != 0) {
message(roi.result$status$msg$message)
stop("No solution")
return(NULL)
}
weights <- roi.result$solution[1:N]
names(weights) <- colnames(R)
out <- list()
out$weights <- weights
out$out <- roi.result$objval
obj_vals <- list()
# Calculate the objective values here so that we can use the moments$mean
# that might be passed in by the user. This will avoid
# the extra call to constrained_objective
port.mean <- -roi.result$objval
names(port.mean) <- "mean"
obj_vals[["mean"]] <- port.mean
out$obj_vals <- obj_vals
# out$call <- call # need to get the call outside of the function
return(out)
}
##### Maximize Return MILP Function #####
#' Maximum Return MILP Optimization
#'
#' This function is called by optimize.portfolio to solve maximum return
#' problems via mixed integer linear programming.
#'
#' @param R xts object of asset returns
#' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
#' @param moments object of moments computed based on objective functions
#' @param target target return value
#' @param solver solver to use
#' @param control list of solver control parameters
#' @author Ross Bennett
maxret_milp_opt <- function(R, constraints, moments, target, solver="glpk", control=NULL){
stopifnot("package:ROI" %in% search() || require("ROI",quietly = TRUE))
plugin <- paste0("ROI.plugin.", solver)
stopifnot(paste0("package:", plugin) %in% search() || require(plugin, quietly=TRUE, character.only=TRUE))
# Check for cleaned returns in moments
if(!is.null(moments$cleanR)) R <- moments$cleanR
# Number of assets
N <- ncol(R)
# Maximum number of positions (non-zero weights)
max_pos <- constraints$max_pos
min_pos <- 1
# Upper and lower bounds on weights
LB <- as.numeric(constraints$min)
UB <- as.numeric(constraints$max)
# Check for target return
if(!is.na(target)){
# We have a target
targetcon <- rbind(c(moments$mean, rep(0, N)),
c(-moments$mean, rep(0, N)))
targetdir <- c("<=", "==")
targetrhs <- c(Inf, -target)
} else {
# No target specified, just maximize
targetcon <- NULL
targetdir <- NULL
targetrhs <- NULL
}
# weight_sum constraint
Amat <- rbind(c(rep(1, N), rep(0, N)),
c(rep(1, N), rep(0, N)))
# Target return constraint
Amat <- rbind(Amat, targetcon)
# Bounds and position limit constraints
Amat <- rbind(Amat, cbind(-diag(N), diag(LB)))
Amat <- rbind(Amat, cbind(diag(N), -diag(UB)))
Amat <- rbind(Amat, c(rep(0, N), rep(-1, N)))
Amat <- rbind(Amat, c(rep(0, N), rep(1, N)))
dir <- c("<=", ">=", targetdir, rep("<=", 2*N), "<=", "<=")
rhs <- c(1, 1, targetrhs, rep(0, 2*N), -min_pos, max_pos)
# Include group constraints
if(try(!is.null(constraints$groups), silent=TRUE)){
n.groups <- length(constraints$groups)
Amat.group <- matrix(0, nrow=n.groups, ncol=N)
k <- 1
l <- 0
for(i in 1:n.groups){
j <- constraints$groups[i]
Amat.group[i, k:(l+j)] <- 1
k <- l + j + 1
l <- k - 1
}
if(is.null(constraints$cLO)) cLO <- rep(-Inf, n.groups)
if(is.null(constraints$cUP)) cUP <- rep(Inf, n.groups)
zeros <- matrix(data=0, nrow=nrow(Amat.group), ncol=ncol(Amat.group))
Amat <- rbind(Amat, cbind(Amat.group, zeros), cbind(-Amat.group, zeros))
dir <- c(dir, rep(">=", (n.groups + n.groups)))
rhs <- c(rhs, constraints$cLO, -constraints$cUP)
}
# Add the factor exposures to Amat, dir, and rhs
if(!is.null(constraints$B)){
t.B <- t(B)
zeros <- matrix(data=0, nrow=nrow(t.B), ncol=ncol(t.B))
Amat <- rbind(Amat, cbind(t.B, zeros), cbind(-t.B, zeros))
dir <- c(dir, rep(">=", 2 * nrow(t.B)))
rhs <- c(rhs, constraints$lower, -constraints$upper)
}
# Only seems to work if I do not specify bounds
# bnds <- V_bound(li=seq.int(1L, 2*m), lb=c(as.numeric(constraints$min), rep(0, m)),
# ui=seq.int(1L, 2*m), ub=c(as.numeric(constraints$max), rep(1, m)))
bnds <- NULL
# Set up the types vector with continuous and binary variables
types <- c(rep("C", N), rep("B", N))
# Set up the linear objective to maximize mean return
ROI_objective <- L_objective(L=c(-moments$mean, rep(0, N)))
# Set up the optimization problem and solve
opt.prob <- OP(objective=ROI_objective,
constraints=L_constraint(L=Amat, dir=dir, rhs=rhs),
bounds=bnds, types=types)
roi.result <- try(ROI_solve(x=opt.prob, solver=solver, control=control), silent=TRUE)
if(inherits(roi.result, "try-error")) stop(paste("No solution found:", roi.result))
# Weights
weights <- roi.result$solution[1:N]
names(weights) <- colnames(R)
# The out object is returned
out <- list()
out$weights <- weights
out$out <- roi.result$objval
obj_vals <- list()
port.mean <- -roi.result$objval
names(port.mean) <- "mean"
obj_vals[["mean"]] <- port.mean
out$obj_vals <- obj_vals
return(out)
}
##### Minimize ETL LP Function #####
#' Minimum ETL LP Optimization
#'
#' This function is called by optimize.portfolio to solve minimum ETL problems.
#'
#' @param R xts object of asset returns
#' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
#' @param moments object of moments computed based on objective functions
#' @param target target return value
#' @param alpha alpha value for ETL/ES/CVaR
#' @param solver solver to use
#' @param control list of solver control parameters
#' @author Ross Bennett
etl_opt <- function(R, constraints, moments, target, alpha, solver="glpk", control=NULL){
stopifnot("package:ROI" %in% search() || require("ROI",quietly = TRUE))
plugin <- paste0("ROI.plugin.", solver)
stopifnot(paste0("package:", plugin) %in% search() || require(plugin, quietly=TRUE, character.only=TRUE))
# Check for cleaned returns in moments
if(!is.null(moments$cleanR)) R <- moments$cleanR
N <- ncol(R)
T <- nrow(R)
# Applying box constraints
LB <- c(as.numeric(constraints$min), rep(0, T), -1)
UB <- c(as.numeric(constraints$max), rep(Inf, T), 1)
bnds <- V_bound(li=seq.int(1L, N+T+1), lb=LB,
ui=seq.int(1L, N+T+1), ub=UB)
# Add this check if mean is not an objective and return is a constraints
if(!is.na(target)){
if(all(moments$mean == 0)){
moments$mean <- colMeans(R)
}
} else {
moments$mean <- rep(0, N)
target <- 0
}
Amat <- cbind(rbind(1, 1, moments$mean, coredata(R)), rbind(0, 0, 0, cbind(diag(T), 1)))
dir.vec <- c(">=","<=",">=",rep(">=",T))
rhs.vec <- c(constraints$min_sum, constraints$max_sum, target ,rep(0, T))
if(try(!is.null(constraints$groups), silent=TRUE)){
n.groups <- length(constraints$groups)
Amat.group <- matrix(0, nrow=n.groups, ncol=N)
for(i in 1:n.groups){
Amat.group[i, constraints$groups[[i]]] <- 1
}
if(is.null(constraints$cLO)) cLO <- rep(-Inf, n.groups)
if(is.null(constraints$cUP)) cUP <- rep(Inf, n.groups)
zeros <- matrix(0, nrow=n.groups, ncol=(T+1))
Amat <- rbind(Amat, cbind(Amat.group, zeros), cbind(-Amat.group, zeros))
dir.vec <- c(dir.vec, rep(">=", (n.groups + n.groups)))
rhs.vec <- c(rhs.vec, constraints$cLO, -constraints$cUP)
}
# Add the factor exposures to Amat, dir, and rhs
if(!is.null(constraints$B)){
t.B <- t(constraints$B)
zeros <- matrix(data=0, nrow=nrow(t.B), ncol=(T+1))
Amat <- rbind(Amat, cbind(t.B, zeros), cbind(-t.B, zeros))
dir.vec <- c(dir.vec, rep(">=", 2 * nrow(t.B)))
rhs.vec <- c(rhs.vec, constraints$lower, -constraints$upper)
}
ROI_objective <- L_objective(c(rep(0,N), rep(1/(alpha*T),T), 1))
opt.prob <- OP(objective=ROI_objective,
constraints=L_constraint(L=Amat, dir=dir.vec, rhs=rhs.vec),
bounds=bnds)
roi.result <- try(ROI_solve(x=opt.prob, solver=solver, control=control), silent=TRUE)
if(inherits(x=roi.result, "try-error")) stop(paste("No solution found:", roi.result))
weights <- roi.result$solution[1:N]
names(weights) <- colnames(R)
out <- list()
out$weights <- weights
out$out <- roi.result$objval
es_names <- c("ES", "ETL", "CVaR")
es_idx <- which(es_names %in% names(moments))
obj_vals <- list()
# Calculate the objective values here so that we can use the moments$mean
# and moments$var that might be passed in by the user. This will avoid
# the extra call to constrained_objective
if(!all(moments$mean == 0)){
port.mean <- as.numeric(sum(weights * moments$mean))
names(port.mean) <- "mean"
obj_vals[["mean"]] <- port.mean
port.es <- roi.result$objval
names(port.es) <- es_names[es_idx]
obj_vals[[es_names[es_idx]]] <- port.es
} else {
port.es <- roi.result$objval
names(port.es) <- es_names[es_idx]
obj_vals[[es_names[es_idx]]] <- port.es
}
out$obj_vals <- obj_vals
return(out)
}
##### Minimize ETL MILP Function #####
#' Minimum ETL MILP Optimization
#'
#' This function is called by optimize.portfolio to solve minimum ETL problems
#' via mixed integer linear programming.
#'
#' @param R xts object of asset returns
#' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
#' @param moments object of moments computed based on objective functions
#' @param target target return value
#' @param alpha alpha value for ETL/ES/CVaR
#' @param solver solver to use
#' @param control list of solver control parameters
#' @author Ross Bennett
etl_milp_opt <- function(R, constraints, moments, target, alpha, solver="glpk", control=NULL){
stopifnot("package:ROI" %in% search() || require("ROI",quietly = TRUE))
plugin <- paste0("ROI.plugin.", solver)
stopifnot(paste0("package:", plugin) %in% search() || require(plugin, quietly=TRUE, character.only=TRUE))
# Check for cleaned returns in moments
if(!is.null(moments$cleanR)) R <- moments$cleanR
# Number of rows
n <- nrow(R)
# Number of columns
m <- ncol(R)
max_sum <- constraints$max_sum
min_sum <- constraints$min_sum
LB <- constraints$min
UB <- constraints$max
max_pos <- constraints$max_pos
min_pos <- 1
moments_mean <- as.numeric(moments$mean)
# A benchmark can be specified in the parma package.
# Leave this in and set to 0 for now
benchmark <- 0
# Check for target return
if(!is.na(target)){
# We have a target
targetcon <- c(moments_mean, rep(0, n+2))
targetdir <- "=="
targetrhs <- target
} else {
# No target specified, just maximize
targetcon <- NULL
targetdir <- NULL
targetrhs <- NULL
}
# Set up initial A matrix
tmpAmat <- cbind(-coredata(R),
matrix(-1, nrow=n, ncol=1),
-diag(n),
matrix(benchmark, nrow=n, ncol=1))
# Add leverage constraints to matrix
tmpAmat <- rbind(tmpAmat, rbind(c(rep(1, m), rep(0, n+2)),
c(rep(1, m), rep(0, n+2))))
# Add target return to matrix
tmpAmat <- rbind(tmpAmat, as.numeric(targetcon))
# This step just adds m rows to the matrix to accept box constraints in the next step
tmpAmat <- cbind(tmpAmat, matrix(0, ncol=m, nrow=dim(tmpAmat)[1]))
# Add lower bound box constraints
tmpAmat <- rbind(tmpAmat, cbind(-diag(m), matrix(0, ncol=n+2, nrow=m), diag(LB)))
# Add upper bound box constraints
tmpAmat <- rbind(tmpAmat, cbind(diag(m), matrix(0, ncol=n+2, nrow=m), diag(-UB)))
# Add rows cardinality constraints
tmpAmat <- rbind(tmpAmat, cbind(matrix(0, ncol=m + n + 2, nrow=1), matrix(-1, ncol=m, nrow=1)))
tmpAmat <- rbind(tmpAmat, cbind(matrix(0, ncol=m + n + 2, nrow=1), matrix(1, ncol=m, nrow=1)))
# Set up the rhs vector
rhs <- c( rep(0, n), min_sum, max_sum, targetrhs, rep(0, 2*m), -min_pos, max_pos)
# Set up the dir vector
dir <- c( rep("<=", n), ">=", "<=", targetdir, rep("<=", 2*m), "<=", "<=")
if(try(!is.null(constraints$groups), silent=TRUE)){
n.groups <- length(constraints$groups)
Amat.group <- matrix(0, nrow=n.groups, ncol=m)
for(i in 1:n.groups){
Amat.group[i, constraints$groups[[i]]] <- 1
}
if(is.null(constraints$cLO)) cLO <- rep(-Inf, n.groups)
if(is.null(constraints$cUP)) cUP <- rep(Inf, n.groups)
zeros <- matrix(0, nrow=n.groups, ncol=(m + n + 2))
tmpAmat <- rbind(tmpAmat, cbind(Amat.group, zeros), cbind(-Amat.group, zeros))
dir <- c(dir, rep(">=", (n.groups + n.groups)))
rhs <- c(rhs, constraints$cLO, -constraints$cUP)
}
# Add the factor exposures to Amat, dir, and rhs
if(!is.null(constraints$B)){
t.B <- t(constraints$B)
zeros <- matrix(data=0, nrow=nrow(t.B), ncol=(m + n + 2))
tmpAmat <- rbind(tmpAmat, cbind(t.B, zeros), cbind(-t.B, zeros))
dir <- c(dir, rep(">=", 2 * nrow(t.B)))
rhs <- c(rhs, constraints$lower, -constraints$upper)
}
# Linear objective vector
ROI_objective <- L_objective(c( rep(0, m), 1, rep(1/n, n) / alpha, 0, rep(0, m)))
# Set up the types vector with continuous and binary variables
types <- c( rep("C", m), "C", rep("C", n), "C", rep("B", m))
bnds <- V_bound( li = 1L:(m + n + 2 + m), lb = c(LB, -1, rep(0, n), 1, rep(0, m)),
ui = 1L:(m + n + 2 + m), ub = c(UB, 1, rep(Inf, n), 1, rep(1, m)))
# Set up the optimization problem and solve
opt.prob <- OP(objective=ROI_objective,
constraints=L_constraint(L=tmpAmat, dir=dir, rhs=rhs),
bounds=bnds, types=types)
roi.result <- try(ROI_solve(x=opt.prob, solver=solver, control=control), silent=TRUE)
if(inherits(roi.result, "try-error")) stop(paste("No solution found:", roi.result))
# The Rglpk solvers status returns an an integer with status information
# about the solution returned: 0 if the optimal solution was found, a
#non-zero value otherwise.
if(roi.result$status$code != 0) {
message("Undefined Solution")
return(NULL)
}
weights <- roi.result$solution[1:m]
names(weights) <- colnames(R)
out <- list()
out$weights <- weights
out$out <- roi.result$objval
es_names <- c("ES", "ETL", "CVaR")
es_idx <- which(es_names %in% names(moments))
obj_vals <- list()
# Calculate the objective values here so that we can use the moments$mean
# and moments$var that might be passed in by the user.
if(!all(moments$mean == 0)){
port.mean <- as.numeric(sum(weights * moments$mean))
names(port.mean) <- "mean"
obj_vals[["mean"]] <- port.mean
port.es <- roi.result$objval
names(port.es) <- es_names[es_idx]
obj_vals[[es_names[es_idx]]] <- port.es
} else {
port.es <- roi.result$objval
names(port.es) <- es_names[es_idx]
obj_vals[[es_names[es_idx]]] <- port.es
}
out$obj_vals <- obj_vals
return(out)
}
##### minimize variance or maximize quadratic utility with turnover constraints #####
#' GMV/QU QP Optimization with Turnover Constraint
#'
#' This function is called by optimize.portfolio to solve minimum variance or
#' maximum quadratic utility problems with turnover constraint
#'
#' @param R xts object of asset returns
#' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
#' @param moments object of moments computed based on objective functions
#' @param lambda risk_aversion parameter
#' @param target target return value
#' @param init_weights initial weights to compute turnover
#' @param solver solver to use
#' @param control list of solver control parameters
#' @author Ross Bennett
gmv_opt_toc <- function(R, constraints, moments, lambda, target, init_weights, solver="quadprog", control=NULL){
# function for minimum variance or max quadratic utility problems
stopifnot("package:corpcor" %in% search() || require("corpcor",quietly = TRUE))
stopifnot("package:ROI" %in% search() || require("ROI", quietly = TRUE))
plugin <- paste0("ROI.plugin.", solver)
stopifnot(paste0("package:", plugin) %in% search() || require(plugin, quietly=TRUE, character.only=TRUE))
# Check for cleaned returns in moments
if(!is.null(moments$cleanR)) R <- moments$cleanR
# Modify the returns matrix. This is done because there are 3 sets of
# variables 1) w.initial, 2) w.buy, and 3) w.sell
R0 <- matrix(0, ncol=ncol(R), nrow=nrow(R))
returns <- cbind(R, R0, R0)
V <- cov(returns)
# number of assets
N <- ncol(R)
# initial weights for solver
if(is.null(init_weights)) init_weights <- rep(1/ N, N)
# check for a target return constraint
if(!is.na(target)) {
# If var is the only objective specified, then moments$mean won't be calculated
if(all(moments$mean==0)){
tmp_means <- colMeans(R)
} else {
tmp_means <- moments$mean
}
} else {
tmp_means <- rep(0, N)
target <- 0
}
Amat <- c(tmp_means, rep(0, 2*N))
dir <- "=="
rhs <- target
meq <- N + 1
# Amat for initial weights
# Amat <- cbind(diag(N), matrix(0, nrow=N, ncol=N*2))
Amat <- rbind(Amat, cbind(diag(N), -1*diag(N), diag(N)))
rhs <- c(rhs, init_weights)
dir <- c(dir, rep("==", N))
# Amat for turnover constraints
Amat <- rbind(Amat, c(rep(0, N), rep(-1, N), rep(-1, N)))
rhs <- c(rhs, -constraints$turnover_target)
dir <- c(dir, ">=")
# Amat for positive weights
Amat <- rbind(Amat, cbind(matrix(0, nrow=N, ncol=N), diag(N), matrix(0, nrow=N, ncol=N)))
rhs <- c(rhs, rep(0, N))
dir <- c(dir, rep(">=", N))
# Amat for negative weights
Amat <- rbind(Amat, cbind(matrix(0, nrow=N, ncol=2*N), diag(N)))
rhs <- c(rhs, rep(0, N))
dir <- c(dir, rep(">=", N))
# Amat for full investment constraint
Amat <- rbind(Amat, rbind(c(rep(1, N), rep(0,2*N)), c(rep(-1, N), rep(0,2*N))))
rhs <- c(rhs, constraints$min_sum, -constraints$max_sum)
dir <- c(dir, ">=", ">=")
# Amat for lower box constraints
Amat <- rbind(Amat, cbind(diag(N), diag(0, N), diag(0, N)))
rhs <- c(rhs, constraints$min)
dir <- c(dir, rep(">=", N))
# Amat for upper box constraints
Amat <- rbind(Amat, cbind(-diag(N), diag(0, N), diag(0, N)))
rhs <- c(rhs, -constraints$max)
dir <- c(dir, rep(">=", N))
# include group constraints
if(try(!is.null(constraints$groups), silent=TRUE)){
n.groups <- length(constraints$groups)
Amat.group <- matrix(0, nrow=n.groups, ncol=N)
zeros <- matrix(0, nrow=n.groups, ncol=N)
for(i in 1:n.groups){
Amat.group[i, constraints$groups[[i]]] <- 1
}
if(is.null(constraints$cLO)) cLO <- rep(-Inf, n.groups)
if(is.null(constraints$cUP)) cUP <- rep(Inf, n.groups)
Amat <- rbind(Amat, cbind(Amat.group, zeros, zeros))
Amat <- rbind(Amat, cbind(-Amat.group, zeros, zeros))
dir <- c(dir, rep(">=", (n.groups + n.groups)))
rhs <- c(rhs, constraints$cLO, -constraints$cUP)
}
# Add the factor exposures to Amat, dir, and rhs
if(!is.null(constraints$B)){
t.B <- t(constraints$B)
zeros <- matrix(0, nrow=nrow(t.B), ncol=ncol(t.B))
Amat <- rbind(Amat, cbind(t.B, zeros, zeros))
Amat <- rbind(Amat, cbind(-t.B, zeros, zeros))
dir <- c(dir, rep(">=", 2 * nrow(t.B)))
rhs <- c(rhs, constraints$lower, -constraints$upper)
}
# Remove the rows of Amat and elements of rhs.vec where rhs is Inf or -Inf
Amat <- Amat[!is.infinite(rhs), ]
rhs <- rhs[!is.infinite(rhs)]
dir <- dir[!is.infinite(rhs)]
ROI_objective <- Q_objective(Q=make.positive.definite(2*lambda*V),
L=rep(-tmp_means, 3))
opt.prob <- OP(objective=ROI_objective,
constraints=L_constraint(L=Amat, dir=dir, rhs=rhs))
roi.result <- try(ROI_solve(x=opt.prob, solver=solver, control=control), silent=TRUE)
if(inherits(roi.result, "try-error")) stop(paste("No solution found:", roi.result))
wts <- roi.result$solution
wts.final <- wts[1:N]
weights <- wts.final
names(weights) <- colnames(R)
out <- list()
out$weights <- weights
out$out <- roi.result$objval
obj_vals <- list()
# Calculate the objective values here so that we can use the moments$mean
# and moments$var that might be passed in by the user.
if(!all(moments$mean == 0)){
port.mean <- as.numeric(sum(weights * moments$mean))
names(port.mean) <- "mean"
obj_vals[["mean"]] <- port.mean
# faster and more efficient way to compute t(w) %*% Sigma %*% w
port.sd <- sqrt(sum(crossprod(weights, moments$var) * weights))
names(port.sd) <- "StdDev"
obj_vals[["StdDev"]] <- port.sd
} else {
# faster and more efficient way to compute t(w) %*% Sigma %*% w
port.sd <- sqrt(sum(crossprod(weights, moments$var) * weights))
names(port.sd) <- "StdDev"
obj_vals[["StdDev"]] <- port.sd
}
out$obj_vals <- obj_vals
return(out)
}
##### minimize variance or maximize quadratic utility with proportional transactioncosts constraints #####
#' GMV/QU QP Optimization with Proportional Transaction Cost Constraint
#'
#' This function is called by optimize.portfolio to solve minimum variance or
#' maximum quadratic utility problems with proportional transaction cost constraint
#'
#' @param R xts object of asset returns
#' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
#' @param moments object of moments computed based on objective functions
#' @param lambda risk_aversion parameter
#' @param target target return value
#' @param init_weights initial weights to compute turnover
#' @param solver solver to use
#' @param control list of solver control parameters
#' @author Ross Bennett
gmv_opt_ptc <- function(R, constraints, moments, lambda, target, init_weights, solver="quadprog", control=NULL){
# function for minimum variance or max quadratic utility problems
# modifying ProportionalCostOpt function from MPO package
stopifnot("package:corpcor" %in% search() || require("corpcor", quietly = TRUE))
stopifnot("package:ROI" %in% search() || require("ROI", quietly = TRUE))
plugin <- paste0("ROI.plugin.", solver)
stopifnot(paste0("package:", plugin) %in% search() || require(plugin, quietly=TRUE, character.only=TRUE))
# Check for cleaned returns in moments
if(!is.null(moments$cleanR)) R <- moments$cleanR
# Modify the returns matrix. This is done because there are 3 sets of
# variables 1) w.initial, 2) w.buy, and 3) w.sell
returns <- cbind(R, R, R)
V <- cov(returns)
# number of assets
N <- ncol(R)
# initial weights for solver
if(is.null(init_weights)) init_weights <- rep(1/ N, N)
# Amat for initial weights
Amat <- cbind(diag(N), matrix(0, nrow=N, ncol=N*2))
rhs <- init_weights
dir <- rep("==", N)
meq <- N
# check for a target return constraint
if(!is.na(target)) {
# If var is the only objective specified, then moments$mean won't be calculated
if(all(moments$mean==0)){
tmp_means <- colMeans(R)
} else {
tmp_means <- moments$mean
}
Amat <- rbind(Amat, rep((1+tmp_means), 3))
dir <- c(dir, "==")
rhs <- c(rhs, (1+target))
meq <- N + 1
}
# Amat for positive weights for w.buy and w.sell
weights.positive <- rbind(matrix(0,ncol=2*N,nrow=N),diag(2*N))
temp.index <- (N*3-N+1):(N*3)
weights.positive[temp.index,] <- -1*weights.positive[temp.index,]
Amat <- rbind(Amat, t(weights.positive))
rhs <- c(rhs, rep(0, 2*N))
# Amat for full investment constraint
ptc <- constraints$ptc
Amat <- rbind(Amat, rbind(c(rep(1, N), (1+ptc), (1-ptc)), -c(rep(1, N), (1+ptc), (1-ptc))))
rhs <- c(rhs, constraints$min_sum, -constraints$max_sum)
dir <- c(dir, ">=", ">=")
# Amat for lower box constraints
Amat <- rbind(Amat, cbind(diag(N), diag(N), diag(N)))
rhs <- c(rhs, constraints$min)
dir <- c(dir, rep(">=", N))
# Amat for upper box constraints
Amat <- rbind(Amat, cbind(-diag(N), -diag(N), -diag(N)))
rhs <- c(rhs, -constraints$max)
dir <- c(dir, rep(">=", N))
# include group constraints
if(try(!is.null(constraints$groups), silent=TRUE)){
n.groups <- length(constraints$groups)
Amat.group <- matrix(0, nrow=n.groups, ncol=N)
for(i in 1:n.groups){
Amat.group[i, constraints$groups[[i]]] <- 1
}
if(is.null(constraints$cLO)) cLO <- rep(-Inf, n.groups)
if(is.null(constraints$cUP)) cUP <- rep(Inf, n.groups)
Amat <- rbind(Amat, cbind(Amat.group, Amat.group, Amat.group))
Amat <- rbind(Amat, cbind(-Amat.group, -Amat.group, -Amat.group))
dir <- c(dir, rep(">=", (n.groups + n.groups)))
rhs <- c(rhs, constraints$cLO, -constraints$cUP)
}
# Add the factor exposures to Amat, dir, and rhs
if(!is.null(constraints$B)){
t.B <- t(constraints$B)
Amat <- rbind(Amat, cbind(t.B, t.B, t.B))
Amat <- rbind(Amat, cbind(-t.B, -t.B, -t.B))
dir <- c(dir, rep(">=", 2 * nrow(t.B)))
rhs <- c(rhs, constraints$lower, -constraints$upper)
}
d <- rep(-moments$mean, 3)
# Remove the rows of Amat and elements of rhs.vec where rhs is Inf or -Inf
Amat <- Amat[!is.infinite(rhs), ]
rhs <- rhs[!is.infinite(rhs)]
dir <- dir[!is.infinite(rhs)]
ROI_objective <- Q_objective(Q=make.positive.definite(2*lambda*V),
L=rep(-moments$mean, 3))
opt.prob <- OP(objective=ROI_objective,
constraints=L_constraint(L=Amat, dir=dir, rhs=rhs))
roi.result <- try(ROI_solve(x=opt.prob, solver=solver, control=control), silent=TRUE)
if(inherits(roi.result, "try-error")) stop(paste("No solution found:", roi.result))
wts <- roi.result$solution
w.buy <- roi.result$solution[(N+1):(2*N)]
w.sell <- roi.result$solution[(2*N+1):(3*N)]
w.total <- init_weights + w.buy + w.sell
wts.final <- wts[(1:N)] + wts[(1+N):(2*N)] + wts[(2*N+1):(3*N)]
weights <- wts.final
names(weights) <- colnames(R)
out <- list()
out$weights <- weights
out$out <- roi.result$objval
obj_vals <- list()
# Calculate the objective values here so that we can use the moments$mean
# and moments$var that might be passed in by the user.
if(!all(moments$mean == 0)){
port.mean <- as.numeric(sum(weights * moments$mean))
names(port.mean) <- "mean"
obj_vals[["mean"]] <- port.mean
# faster and more efficient way to compute t(w) %*% Sigma %*% w
port.sd <- sqrt(sum(crossprod(weights, moments$var) * weights))
names(port.sd) <- "StdDev"
obj_vals[["StdDev"]] <- port.sd
} else {
# faster and more efficient way to compute t(w) %*% Sigma %*% w
port.sd <- sqrt(sum(crossprod(weights, moments$var) * weights))
names(port.sd) <- "StdDev"
obj_vals[["StdDev"]] <- port.sd
}
out$obj_vals <- obj_vals
return(out)
}
# This function uses optimize() to find the target return value that
# results in the maximum starr ratio (mean / ES).
# returns the target return value
mean_etl_opt <- function(R, constraints, moments, alpha, solver, control){
# create a copy of the moments that can be modified
tmp_moments <- moments
# Find the maximum return
if(!is.null(constraints$max_pos)){
max_ret <- maxret_milp_opt(R=R, constraints=constraints, moments=moments, target=NA, solver=solver, control=control)
} else {
max_ret <- maxret_opt(R=R, moments=moments, constraints=constraints, target=NA, solver=solver, control=control)
}
max_mean <- as.numeric(-max_ret$out)
# Find the minimum return
tmp_moments$mean <- -1 * moments$mean
if(!is.null(constraints$max_pos)){
min_ret <- maxret_milp_opt(R=R, constraints=constraints, moments=tmp_moments, target=NA, solver=solver, control=control)
} else {
min_ret <- maxret_opt(R=R, constraints=constraints, moments=tmp_moments, target=NA, solver=solver, control=control)
}
min_mean <- as.numeric(min_ret$out)
# use optimize() to find the target return value that maximizes sharpe ratio
opt <- try(optimize(f=starr_obj_fun, R=R, constraints=constraints,
solver=solver, moments=moments, alpha=alpha,
lower=min_mean, upper=max_mean, control=control,
maximum=TRUE, tol=.Machine$double.eps),
silent=TRUE)
if(inherits(opt, "try-error")){
stop(paste("Objective function failed with message\n", opt))
return(NULL)
}
return(opt$maximum)
}
# Function to calculate the starr ratio.
# Used as the objective function for optimize()
starr_obj_fun <- function(target_return, R, constraints, moments, alpha, solver, control){
if(!is.null(constraints$max_pos)){
opt <- etl_milp_opt(R=R, constraints=constraints, moments=moments,
target=target_return, alpha=alpha, solver=solver,
control=control)
} else {
opt <- etl_opt(R=R, constraints=constraints, moments=moments,
target=target_return, alpha=alpha, solver=solver,
control=control)
}
weights <- matrix(opt$weights, ncol=1)
opt_mean <- sum(weights * moments$mean)
opt_etl <- as.numeric(opt$out)
starr <- opt_mean / opt_etl
return(starr)
}
# This was my implementation of a binary search for the maximum starr ratio
# target return. Better to use optimize() in R rather than my method. -Ross Bennett
# mean_etl_opt <- function(R, constraints, moments, target, alpha, solver="glpk", tol=.Machine$double.eps^0.5, maxit=50){
# # This function returns the target mean return that maximizes mean / etl (i.e. starr)
#
# # if all(moments$mean == 0) then the user did not specify mean as an objective,
# # and we just want to return the target mean return value
# if(all(moments$mean == 0)) return(target)
#
# fmean <- matrix(moments$mean, ncol=1)
#
# # can't use optimize.portfolio here, this function is called inside
# # optimize.portfolio and will throw an error message about nesting too deeply
#
# # Find the maximum return
# if(!is.null(constraints$max_pos)){
# max_ret <- maxret_milp_opt(R=R, constraints=constraints, moments=moments, target=NA, solver=solver)
# } else {
# max_ret <- maxret_opt(R=R, moments=moments, constraints=constraints, target=NA, solver=solver)
# }
# max_mean <- as.numeric(-max_ret$out)
#
# # Find the starr at the maximum etl portfolio
# if(!is.null(constraints$max_pos)){
# ub_etl <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=max_mean, alpha=alpha, solver=solver)
# } else {
# ub_etl <- etl_opt(R=R, constraints=constraints, moments=moments, target=max_mean, alpha=alpha, solver=solver)
# }
# ub_weights <- matrix(ub_etl$weights, ncol=1)
# ub_mean <- as.numeric(t(ub_weights) %*% fmean)
# ub_etl <- as.numeric(ub_etl$out)
# # starr at the upper bound
# ub_starr <- ub_mean / ub_etl
# if(is.infinite(ub_starr)) stop("Inf value for STARR, objective value is 0")
#
# # cat("ub_mean", ub_mean, "\n")
# # cat("ub_etl", ub_etl, "\n")
# # cat("ub_starr", ub_starr, "\n")
#
# # Find the starr at the minimum etl portfolio
# if(!is.null(constraints$max_pos)){
# lb_etl <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=NA, alpha=alpha, solver=solver)
# } else {
# lb_etl <- etl_opt(R=R, constraints=constraints, moments=moments, target=NA, alpha=alpha, solver=solver)
# }
# lb_weights <- matrix(lb_etl$weights)
# lb_mean <- as.numeric(t(lb_weights) %*% fmean)
# lb_etl <- as.numeric(lb_etl$out)
#
# # starr at the lower bound
# lb_starr <- lb_mean / lb_etl
# # if(is.infinite(lb_starr)) stop("Inf value for STARR, objective value is 0")
#
# # set lb_starr equal to 0, should this be a negative number like -1e6?
# # the lb_* values will be 0 for a dollar-neutral strategy so we need to reset the values
# if(is.na(lb_starr) | is.infinite(lb_starr)) lb_starr <- 0
#
# # cat("lb_mean", lb_mean, "\n")
# # cat("lb_etl", lb_etl, "\n")
# # cat("lb_starr", lb_starr, "\n")
#
# # want to find the return that maximizes mean / etl
# i <- 1
# while((abs(ub_starr - lb_starr) > tol) & (i < maxit)){
# # bisection method to find the maximum mean / etl
#
# # print(i)
# # cat("ub_starr", ub_starr, "\n")
# # cat("lb_starr", lb_starr, "\n")
# # print("**********")
# # Find the starr at the mean return midpoint
# new_ret <- (lb_mean + ub_mean) / 2
# if(!is.null(constraints$max_pos)){
# mid <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha, solver=solver)
# } else {
# mid <- etl_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha, solver=solver)
# }
# # print(mid)
# mid_weights <- matrix(mid$weights, ncol=1)
# mid_mean <- as.numeric(t(mid_weights) %*% fmean)
# mid_etl <- as.numeric(mid$out)
# mid_starr <- mid_mean / mid_etl
# # the mid_* values MIGHT be 0 for a dollar-neutral strategy so we need to reset the values
# # if(is.na(mid_starr) | is.infinite(mid_starr)) mid_starr <- 0
# # tmp_starr <- mid_starr
#
# # cat("mid_mean", mid_mean, "\n")
# # cat("mid_etl", mid_etl, "\n")
# # cat("mid_starr", mid_starr, "\n")
#
# if(mid_starr > ub_starr){
# # if mid_starr > ub_starr then mid_starr becomes the new upper bound
# ub_mean <- mid_mean
# ub_starr <- mid_starr
# new_ret <- (lb_mean + ub_mean) / 2
# if(!is.null(constraints$max_pos)){
# mid <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha, solver=solver)
# } else {
# mid <- etl_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha, solver=solver)
# }
# mid_weights <- matrix(mid$weights, ncol=1)
# mid_mean <- as.numeric(t(mid_weights) %*% fmean)
# mid_etl <- as.numeric(mid$out)
# mid_starr <- mid_mean / mid_etl
# # the mid_* values MIGHT be 0 for a dollar-neutral strategy so we need to reset the values
# # if(is.na(mid_starr) | is.infinite(mid_starr)) mid_starr <- 0
# } else if(mid_starr > lb_starr){
# # if mid_starr > lb_starr then mid_starr becomes the new lower bound
# lb_mean <- mid_mean
# lb_starr <- mid_starr
# new_ret <- (lb_mean + ub_mean) / 2
# if(!is.null(constraints$max_pos)){
# mid <- etl_milp_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha, solver=solver)
# } else {
# mid <- etl_opt(R=R, constraints=constraints, moments=moments, target=new_ret, alpha=alpha, solver=solver)
# }
# mid_weights <- matrix(mid$weights, ncol=1)
# mid_mean <- as.numeric(t(mid_weights) %*% fmean)
# mid_etl <- as.numeric(mid$out)
# mid_starr <- mid_mean / mid_etl
# # the mid_* values MIGHT be 0 for a dollar-neutral strategy so we need to reset the values
# # if(is.na(mid_starr) | is.infinite(mid_starr)) mid_starr <- 0
# }
# i <- i + 1
# }
# return(new_ret)
# }
# Function to calculate the sharpe ratio.
# Used as the objective function for optimize()
sharpe_obj_fun <- function(target_return, R, constraints, moments, lambda_hhi=NULL, conc_groups=NULL, solver="quadprog", control=control){
opt <- gmv_opt(R=R, constraints=constraints, moments=moments, lambda=1,
target=target_return, lambda_hhi=lambda_hhi,
conc_groups=conc_groups, solver=solver, control=control)
weights <- opt$weights
# opt_mean <- as.numeric(t(weights) %*% matrix(moments$mean, ncol=1))
opt_mean <- sum(weights * moments$mean)
# opt_sd <- as.numeric(sqrt(t(weights) %*% moments$var %*% weights))
# faster and more efficient way to compute t(w) %*% Sigma %*% w
opt_sd <- sqrt(sum(crossprod(weights, moments$var) * weights))
opt_sr <- opt_mean / opt_sd
return(opt_sr)
}
# This function uses optimize() to find the target return value that
# results in the maximum sharpe ratio (mean / sd).
# returns the target return value
max_sr_opt <- function(R, constraints, moments, lambda_hhi, conc_groups, solver, control){
# create a copy of the moments that can be modified
tmp_moments <- moments
# Find the maximum return
max_ret <- maxret_opt(R=R, moments=moments, constraints=constraints,
target=NA, solver=solver, control=control)
max_mean <- as.numeric(-max_ret$out)
# Find the minimum return
tmp_moments$mean <- -1 * moments$mean
min_ret <- maxret_opt(R=R, moments=tmp_moments, constraints=constraints,
target=NA, solver=solver, control=control)
min_mean <- as.numeric(min_ret$out)
# use optimize() to find the target return value that maximizes sharpe ratio
opt <- try(optimize(f=sharpe_obj_fun, R=R, constraints=constraints,
solver=solver, lambda_hhi=lambda_hhi,
conc_groups=conc_groups, moments=moments, control=control,
lower=min_mean, upper=max_mean,
maximum=TRUE, tol=.Machine$double.eps),
silent=TRUE)
if(inherits(opt, "try-error")){
stop(paste("Objective function failed with message\n", opt))
return(NULL)
}
return(opt$maximum)
}
# This was my implementation of a binary search for the maximum sharpe ratio
# target return. Better to use optimize() in R rather than my method. -Ross Bennett
# max_sr_opt <- function(R, constraints, moments, lambda, target, lambda_hhi, conc_groups, solver="quadprog", tol=.Machine$double.eps^0.5, maxit=50){
# # This function returns the target mean return that maximizes mean / sd (i.e. sharpe ratio)
#
# # get the forecast mean from moments
# fmean <- matrix(moments$mean, ncol=1)
#
# # Find the maximum return
# max_ret <- maxret_opt(R=R, moments=moments, constraints=constraints, target=NA)
# max_mean <- as.numeric(-max_ret$out)
#
# # Calculate the sr at the maximum mean return portfolio
# ub_weights <- matrix(max_ret$weights, ncol=1)
# ub_mean <- max_mean
# ub_sd <- as.numeric(sqrt(t(ub_weights) %*% moments$var %*% ub_weights))
# # sr at the upper bound
# ub_sr <- ub_mean / ub_sd
#
# # Calculate the sr at the miminum var portfolio
# tmpmoments <- moments
# tmpmoments$mean <- rep(0, length(moments$mean))
# lb_sr <- gmv_opt(R=R, constraints=constraints, moments=tmpmoments, lambda=1, target=NA, lambda_hhi=lambda_hhi, conc_groups=conc_groups, solver=solver)
# lb_weights <- matrix(lb_sr$weights)
# lb_mean <- as.numeric(t(lb_weights) %*% fmean)
# lb_sd <- as.numeric(sqrt(t(lb_weights) %*% moments$var %*% lb_weights))
# # sr at the lower bound
# lb_sr <- lb_mean / lb_sd
#
# # cat("lb_mean:", lb_mean, "\n")
# # cat("ub_mean:", ub_mean, "\n")
# # print("**********")
#
# # want to find the return that maximizes mean / sd
# i <- 1
# while((abs(ub_sr - lb_sr) > tol) & (i < maxit)){
# # bisection method to find the maximum mean / sd
#
# # Find the starr at the mean return midpoint
# new_ret <- (lb_mean + ub_mean) / 2
# mid <- gmv_opt(R=R, constraints=constraints, moments=tmpmoments, lambda=1, target=new_ret, lambda_hhi=lambda_hhi, conc_groups=conc_groups, solver=solver)
# mid_weights <- matrix(mid$weights, ncol=1)
# mid_mean <- as.numeric(t(mid_weights) %*% fmean)
# mid_sd <- as.numeric(sqrt(t(mid_weights) %*% moments$var %*% mid_weights))
# mid_sr <- mid_mean / mid_sd
# # tmp_sr <- mid_sr
#
# # print(i)
# # cat("new_ret:", new_ret, "\n")
# # cat("mid_sr:", mid_sr, "\n")
# # print("**********")
#
# if(mid_sr > ub_sr){
# # if mid_sr > ub_sr then mid_sr becomes the new upper bound
# ub_mean <- mid_mean
# ub_sr <- mid_sr
# new_ret <- (lb_mean + ub_mean) / 2
# mid <- gmv_opt(R=R, constraints=constraints, moments=tmpmoments, lambda=1, target=new_ret, lambda_hhi=lambda_hhi, conc_groups=conc_groups, solver=solver)
# mid_weights <- matrix(mid$weights, ncol=1)
# mid_mean <- as.numeric(t(mid_weights) %*% fmean)
# mid_sd <- as.numeric(sqrt(t(mid_weights) %*% moments$var %*% mid_weights))
# mid_sr <- mid_mean / mid_sd
# } else if(mid_sr > lb_sr){
# # if mid_sr > lb_sr then mid_sr becomes the new lower bound
# lb_mean <- mid_mean
# lb_sr <- mid_sr
# new_ret <- (lb_mean + ub_mean) / 2
# mid <- gmv_opt(R=R, constraints=constraints, moments=tmpmoments, lambda=1, target=new_ret, lambda_hhi=lambda_hhi, conc_groups=conc_groups, solver=solver)
# mid_weights <- matrix(mid$weights, ncol=1)
# mid_mean <- as.numeric(t(mid_weights) %*% fmean)
# mid_sd <- as.numeric(sqrt(t(mid_weights) %*% moments$var %*% mid_weights))
# mid_sr <- mid_mean / mid_sd
# }
# i <- i + 1
# }
# return(new_ret)
# }
###############################################################################
# R (http://r-project.org/) Numeric Methods for Optimization of Portfolios
#
# Copyright (c) 2004-2014 Brian G. Peterson, Peter Carl, Ross Bennett, Kris Boudt
#
# This library is distributed under the terms of the GNU Public License (GPL)
# for full details see the file COPYING
#
# $Id$
#
###############################################################################
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.