[Returnanalytics-commits] r3489 - in pkg/PortfolioAnalytics: R sandbox
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jul 30 00:55:03 CEST 2014
Author: rossbennett34
Date: 2014-07-30 00:55:01 +0200 (Wed, 30 Jul 2014)
New Revision: 3489
Added:
pkg/PortfolioAnalytics/sandbox/leverageQP.R
Modified:
pkg/PortfolioAnalytics/R/optFUN.R
pkg/PortfolioAnalytics/R/optimize.portfolio.R
Log:
adding qp formulation for leverage constraint
Modified: pkg/PortfolioAnalytics/R/optFUN.R
===================================================================
--- pkg/PortfolioAnalytics/R/optFUN.R 2014-07-29 19:25:09 UTC (rev 3488)
+++ pkg/PortfolioAnalytics/R/optFUN.R 2014-07-29 22:55:01 UTC (rev 3489)
@@ -965,6 +965,162 @@
return(out)
}
+##### minimize variance or maximize quadratic utility with leverage 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 a leverage 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 solver solver to use
+#' @param control list of solver control parameters
+#' @author Ross Bennett
+gmv_opt_leverage <- function(R, constraints, moments, lambda, target, 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)
+
+ # 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
+
+ # separate the weights into w, w+, and w-
+ # w - w+ + w- = 0
+ Amat <- rbind(Amat, cbind(diag(N), -1*diag(N), diag(N)))
+ rhs <- c(rhs, rep(0, N))
+ dir <- c(dir, rep("==", N))
+
+ # Amat for leverage constraints
+ Amat <- rbind(Amat, c(rep(0, N), rep(-1, N), rep(-1, N)))
+ rhs <- c(rhs, -constraints$leverage)
+ 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)
+}
+
# This function uses optimize() to find the target return value that
# results in the maximum starr ratio (mean / ES).
# returns the target return value
Modified: pkg/PortfolioAnalytics/R/optimize.portfolio.R
===================================================================
--- pkg/PortfolioAnalytics/R/optimize.portfolio.R 2014-07-29 19:25:09 UTC (rev 3488)
+++ pkg/PortfolioAnalytics/R/optimize.portfolio.R 2014-07-29 22:55:01 UTC (rev 3489)
@@ -887,11 +887,12 @@
}
# Minimize variance if the only objective specified is variance
# Maximize Quadratic Utility if var and mean are specified as objectives
- if(!is.null(constraints$turnover_target) | !is.null(constraints$ptc)){
+ if(!is.null(constraints$turnover_target) | !is.null(constraints$ptc) | !is.null(constraints$leverage)){
if(!is.null(constraints$turnover_target) & !is.null(constraints$ptc)){
warning("Turnover and proportional transaction cost constraints detected, only running optimization for turnover constraint.")
constraints$ptc <- NULL
}
+ # turnover constraint
if(!is.null(constraints$turnover_target) & is.null(constraints$ptc)){
qp_result <- gmv_opt_toc(R=R, constraints=constraints, moments=moments, lambda=lambda, target=target, init_weights=portfolio$assets, solver=solver, control=control)
weights <- qp_result$weights
@@ -899,6 +900,7 @@
obj_vals <- qp_result$obj_vals
out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=qp_result$out, call=call)
}
+ # proportional transaction costs constraint
if(!is.null(constraints$ptc) & is.null(constraints$turnover_target)){
qp_result <- gmv_opt_ptc(R=R, constraints=constraints, moments=moments, lambda=lambda, target=target, init_weights=portfolio$assets, solver=solver, control=control)
weights <- qp_result$weights
@@ -906,6 +908,14 @@
obj_vals <- qp_result$obj_vals
out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=qp_result$out, call=call)
}
+ # leverage constraint
+ if(!is.null(constraints$leverage)){
+ qp_result <- gmv_opt_leverage(R=R, constraints=constraints, moments=moments, lambda=lambda, target=target, solver=solver, control=control)
+ weights <- qp_result$weights
+ # obj_vals <- constrained_objective(w=weights, R=R, portfolio, trace=TRUE, normalize=FALSE)$objective_measures
+ obj_vals <- qp_result$obj_vals
+ out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=qp_result$out, call=call)
+ }
} else {
# if(hasArg(ef)) ef=match.call(expand.dots=TRUE)$ef else ef=FALSE
if(hasArg(maxSR)) maxSR=match.call(expand.dots=TRUE)$maxSR else maxSR=FALSE
Added: pkg/PortfolioAnalytics/sandbox/leverageQP.R
===================================================================
--- pkg/PortfolioAnalytics/sandbox/leverageQP.R (rev 0)
+++ pkg/PortfolioAnalytics/sandbox/leverageQP.R 2014-07-29 22:55:01 UTC (rev 3489)
@@ -0,0 +1,83 @@
+
+# leverage constrained minimum variance QP
+
+library(PortfolioAnalytics)
+library(corpcor)
+library(quadprog)
+
+data(edhec)
+R <- edhec[, 1:10]
+
+N <- ncol(R)
+leverage <- 1.6
+min_sum <- 0.99
+max_sum <- 1.01
+min_box <- rep(-0.3, N)
+max_box <- rep(1, N)
+
+R0 <- matrix(0, ncol=ncol(R), nrow=nrow(R))
+returns <- cbind(R, R0, R0)
+V <- corpcor::make.positive.definite(cov(returns))
+
+# separate the weights into w, w+, and w-
+# w - w+ + w- = 0
+Amat <- cbind(diag(N), -diag(N), diag(N))
+rhs <- rep(0, N)
+
+# leverage constraint
+# w+ + w- <= leverage
+Amat <- rbind(Amat, c(rep(0, N), rep(-1, N), rep(-1, N)))
+rhs <- c(rhs, -leverage)
+
+# w+ >= 0
+Amat <- rbind(Amat, cbind(diag(0, N), diag(N), diag(0, N)))
+rhs <- c(rhs, rep(0, N))
+
+# w- >= 0
+Amat <- rbind(Amat, cbind(diag(0, N), diag(0, N), diag(N)))
+rhs <- c(rhs, rep(0, N))
+
+# w^T 1 >= min_sum
+Amat <- rbind(Amat, c(rep(1, N), rep(0, N), rep(0, N)))
+rhs <- c(rhs, min_sum)
+
+# w^T 1 <= max_sum
+Amat <- rbind(Amat, c(rep(-1, N), rep(0, N), rep(0, N)))
+rhs <- c(rhs, -max_sum)
+
+# lower box constraints
+Amat <- rbind(Amat, cbind(diag(N), diag(0, N), diag(0, N)))
+rhs <- c(rhs, min_box)
+
+# upper box constraints
+Amat <- rbind(Amat, cbind(-diag(N), diag(0, N), diag(0, N)))
+rhs <- c(rhs, -max_box)
+
+sol <- solve.QP(Dmat=V, dvec=rep(0, 3*N), Amat=t(Amat), bvec=rhs, meq=N)
+sol
+
+weights <- sol$solution[1:N]
+weights
+sum(weights)
+sum(abs(weights)) <= leverage
+
+
+#' This script demonstrates how to solve a constrained portfolio optimization
+#' problem to minimize standard deviation.
+
+#' Load the package and data
+library(PortfolioAnalytics)
+data(edhec)
+R <- edhec[, 1:10]
+funds <- colnames(R)
+
+#' Construct initial portfolio with basic constraints.
+init.portf <- portfolio.spec(assets=funds)
+init.portf <- add.constraint(portfolio=init.portf, type="full_investment")
+init.portf <- add.constraint(portfolio=init.portf, type="box", min=-0.3, max=1)
+init.portf <- add.constraint(portfolio=init.portf, type="leverage_exposure", leverage=1.6)
+init.portf <- add.objective(portfolio=init.portf, type="risk", name="StdDev")
+
+opt <- optimize.portfolio(R, init.portf, optimize_method="ROI")
+opt$weights
+
More information about the Returnanalytics-commits
mailing list