[Returnanalytics-commits] r2899 - in pkg/PortfolioAnalytics: R man sandbox
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Aug 27 04:53:19 CEST 2013
Author: rossbennett34
Date: 2013-08-27 04:53:19 +0200 (Tue, 27 Aug 2013)
New Revision: 2899
Added:
pkg/PortfolioAnalytics/man/etl_milp_opt.Rd
pkg/PortfolioAnalytics/man/etl_opt.Rd
pkg/PortfolioAnalytics/man/gmv_opt.Rd
pkg/PortfolioAnalytics/man/gmv_opt_toc.Rd
pkg/PortfolioAnalytics/man/maxret_milp_opt.Rd
pkg/PortfolioAnalytics/man/maxret_opt.Rd
Modified:
pkg/PortfolioAnalytics/R/optFUN.R
pkg/PortfolioAnalytics/R/optimize.portfolio.R
pkg/PortfolioAnalytics/sandbox/testing_turnover.gmv.R
Log:
Adding function for gmv/qu with turnover constraints. Adding documentation for optimization QP/LP optimization functions.
Modified: pkg/PortfolioAnalytics/R/optFUN.R
===================================================================
--- pkg/PortfolioAnalytics/R/optFUN.R 2013-08-26 19:17:35 UTC (rev 2898)
+++ pkg/PortfolioAnalytics/R/optFUN.R 2013-08-27 02:53:19 UTC (rev 2899)
@@ -1,5 +1,15 @@
##### GMV and QU QP Function #####
+#' Optimization function to solve minimum variance or maximum quadratic utility problems
+#'
+#' This function is called by optimize.portfolio to solve minimum variance or maximum quadratic utility problems
+#'
+#' @param R xts object of asset returns
+#' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
+#' @param moments object of moments computed based on objective functions
+#' @param lambda risk_aversion parameter
+#' @param target target return value
+#' @author Ross Bennett
gmv_opt <- function(R, constraints, moments, lambda, target){
N <- ncol(R)
@@ -66,6 +76,15 @@
}
##### Maximize Return LP Function #####
+#' Optimization function to solve minimum variance or maximum quadratic utility problems
+#'
+#' This function is called by optimize.portfolio to solve minimum variance or maximum quadratic utility problems
+#'
+#' @param R xts object of asset returns
+#' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
+#' @param moments object of moments computed based on objective functions
+#' @param target target return value
+#' @author Ross Bennett
maxret_opt <- function(R, moments, constraints, target){
N <- ncol(R)
@@ -137,6 +156,15 @@
}
##### Maximize Return MILP Function #####
+#' Optimization function to solve minimum variance or maximum quadratic utility problems
+#'
+#' This function is called by optimize.portfolio to solve minimum variance or maximum quadratic utility problems
+#'
+#' @param R xts object of asset returns
+#' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
+#' @param moments object of moments computed based on objective functions
+#' @param target target return value
+#' @author Ross Bennett
maxret_milp_opt <- function(R, constraints, moments, target){
N <- ncol(R)
@@ -226,6 +254,16 @@
}
##### Minimize ETL LP Function #####
+#' Optimization function to solve minimum variance or maximum quadratic utility problems
+#'
+#' This function is called by optimize.portfolio to solve minimum variance or maximum quadratic utility problems
+#'
+#' @param R xts object of asset returns
+#' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
+#' @param moments object of moments computed based on objective functions
+#' @param target target return value
+#' @param alpha alpha value for ETL/ES/CVaR
+#' @author Ross Bennett
etl_opt <- function(R, constraints, moments, target, alpha){
N <- ncol(R)
@@ -277,6 +315,16 @@
}
##### Minimize ETL MILP Function #####
+#' Optimization function to solve minimum variance or maximum quadratic utility problems
+#'
+#' This function is called by optimize.portfolio to solve minimum variance or maximum quadratic utility problems
+#'
+#' @param R xts object of asset returns
+#' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
+#' @param moments object of moments computed based on objective functions
+#' @param target target return value
+#' @param alpha alpha value for ETL/ES/CVaR
+#' @author Ross Bennett
etl_milp_opt <- function(R, constraints, moments, target, alpha){
# Number of rows
@@ -389,3 +437,132 @@
#out$call <- call # add this outside of here, this function doesn't have the call
return(out)
}
+
+##### minimize variance or maximize quadratic utility with turnover constraints #####
+#' Optimization function to solve minimum variance or maximum quadratic utility problems
+#'
+#' This function is called by optimize.portfolio to solve minimum variance or maximum quadratic utility problems
+#'
+#' @param R xts object of asset returns
+#' @param constraints object of constraints in the portfolio object extracted with \code{get_constraints}
+#' @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
+#' @author Ross Bennett
+gmv_opt_toc <- function(R, constraints, moments, lambda, target, init_weights){
+ # function for minimum variance or max quadratic utility problems
+
+ # 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 <- 4
+
+ # 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(tmp_means, 3))
+ dir <- c(dir, "==")
+ rhs <- c(rhs, target)
+ meq <- 5
+ }
+
+ # Amat for full investment constraint
+ Amat <- rbind(Amat, rbind(rep(1, N*3), rep(-1, N*3)))
+ 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))
+
+ # Amat for turnover constraints
+ Amat <- rbind(Amat, c(rep(0, N), rep(-1, N), rep(1, N)))
+ rhs <- c(rhs, -constraints$toc)
+ 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))
+
+ # 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)
+
+ qp.result <- try(solve.QP(Dmat=make.positive.definite(2*lambda*V),
+ dvec=d, Amat=t(Amat), bvec=rhs, meq=meq), silent=TRUE)
+ if(inherits(qp.result, "try-error")) stop("No solution found, consider adjusting constraints.")
+
+ wts <- qp.result$solution
+ 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 <- qp.result$val
+ return(out)
+
+ # TODO
+ # Get this working with ROI
+
+ # Not getting solution using ROI
+ # set up the quadratic objective
+ # 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 <- ROI_solve(x=opt.prob, solver="quadprog")
+}
+
Modified: pkg/PortfolioAnalytics/R/optimize.portfolio.R
===================================================================
--- pkg/PortfolioAnalytics/R/optimize.portfolio.R 2013-08-26 19:17:35 UTC (rev 2898)
+++ pkg/PortfolioAnalytics/R/optimize.portfolio.R 2013-08-27 02:53:19 UTC (rev 2899)
@@ -703,10 +703,17 @@
if("var" %in% names(moments)){
# Minimize variance if the only objective specified is variance
# Maximize Quadratic Utility if var and mean are specified as objectives
- roi_result <- gmv_opt(R=R, constraints=constraints, moments=moments, lambda=lambda, target=target)
- weights <- roi_result$weights
- obj_vals <- constrained_objective(w=weights, R=R, portfolio, trace=TRUE, normalize=FALSE)$objective_measures
- out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=roi_result$out, call=call)
+ if(!is.null(constraints$turnover_target)){
+ qp_result <- gmv_opt_toc(R=R, constraints=constraints, moments=moments, lambda=lambda, target=target, init_weights=portfolio$assets)
+ weights <- qp_result$weights
+ obj_vals <- constrained_objective(w=weights, R=R, portfolio, trace=TRUE, normalize=FALSE)$objective_measures
+ out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=roi_result$out, call=call)
+ } else {
+ roi_result <- gmv_opt(R=R, constraints=constraints, moments=moments, lambda=lambda, target=target)
+ weights <- roi_result$weights
+ obj_vals <- constrained_objective(w=weights, R=R, portfolio, trace=TRUE, normalize=FALSE)$objective_measures
+ out <- list(weights=weights, objective_measures=obj_vals, opt_values=obj_vals, out=roi_result$out, call=call)
+ }
}
if(length(names(moments)) == 1 & "mean" %in% names(moments)) {
# Maximize return if the only objective specified is mean
Added: pkg/PortfolioAnalytics/man/etl_milp_opt.Rd
===================================================================
--- pkg/PortfolioAnalytics/man/etl_milp_opt.Rd (rev 0)
+++ pkg/PortfolioAnalytics/man/etl_milp_opt.Rd 2013-08-27 02:53:19 UTC (rev 2899)
@@ -0,0 +1,27 @@
+\name{etl_milp_opt}
+\alias{etl_milp_opt}
+\title{Optimization function to solve minimum variance or maximum quadratic utility problems}
+\usage{
+ etl_milp_opt(R, constraints, moments, target, alpha)
+}
+\arguments{
+ \item{R}{xts object of asset returns}
+
+ \item{constraints}{object of constraints in the portfolio
+ object extracted with \code{get_constraints}}
+
+ \item{moments}{object of moments computed based on
+ objective functions}
+
+ \item{target}{target return value}
+
+ \item{alpha}{alpha value for ETL/ES/CVaR}
+}
+\description{
+ This function is called by optimize.portfolio to solve
+ minimum variance or maximum quadratic utility problems
+}
+\author{
+ Ross Bennett
+}
+
Added: pkg/PortfolioAnalytics/man/etl_opt.Rd
===================================================================
--- pkg/PortfolioAnalytics/man/etl_opt.Rd (rev 0)
+++ pkg/PortfolioAnalytics/man/etl_opt.Rd 2013-08-27 02:53:19 UTC (rev 2899)
@@ -0,0 +1,27 @@
+\name{etl_opt}
+\alias{etl_opt}
+\title{Optimization function to solve minimum variance or maximum quadratic utility problems}
+\usage{
+ etl_opt(R, constraints, moments, target, alpha)
+}
+\arguments{
+ \item{R}{xts object of asset returns}
+
+ \item{constraints}{object of constraints in the portfolio
+ object extracted with \code{get_constraints}}
+
+ \item{moments}{object of moments computed based on
+ objective functions}
+
+ \item{target}{target return value}
+
+ \item{alpha}{alpha value for ETL/ES/CVaR}
+}
+\description{
+ This function is called by optimize.portfolio to solve
+ minimum variance or maximum quadratic utility problems
+}
+\author{
+ Ross Bennett
+}
+
Added: pkg/PortfolioAnalytics/man/gmv_opt.Rd
===================================================================
--- pkg/PortfolioAnalytics/man/gmv_opt.Rd (rev 0)
+++ pkg/PortfolioAnalytics/man/gmv_opt.Rd 2013-08-27 02:53:19 UTC (rev 2899)
@@ -0,0 +1,27 @@
+\name{gmv_opt}
+\alias{gmv_opt}
+\title{Optimization function to solve minimum variance or maximum quadratic utility problems}
+\usage{
+ gmv_opt(R, constraints, moments, lambda, target)
+}
+\arguments{
+ \item{R}{xts object of asset returns}
+
+ \item{constraints}{object of constraints in the portfolio
+ object extracted with \code{get_constraints}}
+
+ \item{moments}{object of moments computed based on
+ objective functions}
+
+ \item{lambda}{risk_aversion parameter}
+
+ \item{target}{target return value}
+}
+\description{
+ This function is called by optimize.portfolio to solve
+ minimum variance or maximum quadratic utility problems
+}
+\author{
+ Ross Bennett
+}
+
Added: pkg/PortfolioAnalytics/man/gmv_opt_toc.Rd
===================================================================
--- pkg/PortfolioAnalytics/man/gmv_opt_toc.Rd (rev 0)
+++ pkg/PortfolioAnalytics/man/gmv_opt_toc.Rd 2013-08-27 02:53:19 UTC (rev 2899)
@@ -0,0 +1,30 @@
+\name{gmv_opt_toc}
+\alias{gmv_opt_toc}
+\title{Optimization function to solve minimum variance or maximum quadratic utility problems}
+\usage{
+ gmv_opt_toc(R, constraints, moments, lambda, target,
+ init_weights)
+}
+\arguments{
+ \item{R}{xts object of asset returns}
+
+ \item{constraints}{object of constraints in the portfolio
+ object extracted with \code{get_constraints}}
+
+ \item{moments}{object of moments computed based on
+ objective functions}
+
+ \item{lambda}{risk_aversion parameter}
+
+ \item{target}{target return value}
+
+ \item{init_weights}{initial weights to compute turnover}
+}
+\description{
+ This function is called by optimize.portfolio to solve
+ minimum variance or maximum quadratic utility problems
+}
+\author{
+ Ross Bennett
+}
+
Added: pkg/PortfolioAnalytics/man/maxret_milp_opt.Rd
===================================================================
--- pkg/PortfolioAnalytics/man/maxret_milp_opt.Rd (rev 0)
+++ pkg/PortfolioAnalytics/man/maxret_milp_opt.Rd 2013-08-27 02:53:19 UTC (rev 2899)
@@ -0,0 +1,25 @@
+\name{maxret_milp_opt}
+\alias{maxret_milp_opt}
+\title{Optimization function to solve minimum variance or maximum quadratic utility problems}
+\usage{
+ maxret_milp_opt(R, constraints, moments, target)
+}
+\arguments{
+ \item{R}{xts object of asset returns}
+
+ \item{constraints}{object of constraints in the portfolio
+ object extracted with \code{get_constraints}}
+
+ \item{moments}{object of moments computed based on
+ objective functions}
+
+ \item{target}{target return value}
+}
+\description{
+ This function is called by optimize.portfolio to solve
+ minimum variance or maximum quadratic utility problems
+}
+\author{
+ Ross Bennett
+}
+
Added: pkg/PortfolioAnalytics/man/maxret_opt.Rd
===================================================================
--- pkg/PortfolioAnalytics/man/maxret_opt.Rd (rev 0)
+++ pkg/PortfolioAnalytics/man/maxret_opt.Rd 2013-08-27 02:53:19 UTC (rev 2899)
@@ -0,0 +1,25 @@
+\name{maxret_opt}
+\alias{maxret_opt}
+\title{Optimization function to solve minimum variance or maximum quadratic utility problems}
+\usage{
+ maxret_opt(R, moments, constraints, target)
+}
+\arguments{
+ \item{R}{xts object of asset returns}
+
+ \item{constraints}{object of constraints in the portfolio
+ object extracted with \code{get_constraints}}
+
+ \item{moments}{object of moments computed based on
+ objective functions}
+
+ \item{target}{target return value}
+}
+\description{
+ This function is called by optimize.portfolio to solve
+ minimum variance or maximum quadratic utility problems
+}
+\author{
+ Ross Bennett
+}
+
Modified: pkg/PortfolioAnalytics/sandbox/testing_turnover.gmv.R
===================================================================
--- pkg/PortfolioAnalytics/sandbox/testing_turnover.gmv.R 2013-08-26 19:17:35 UTC (rev 2898)
+++ pkg/PortfolioAnalytics/sandbox/testing_turnover.gmv.R 2013-08-27 02:53:19 UTC (rev 2899)
@@ -2,132 +2,156 @@
# script to solve the gmv optimization problem with turnover constraints using quadprog or ROI
library(PortfolioAnalytics)
-library(PerformanceAnalytics)
+library(ROI)
+library(ROI.plugin.quadprog)
library(quadprog)
library(corpcor)
-# TODO Add documentation for function
+data(edhec)
+R <- edhec[, 1:4]
+init.weights <- rep(1/4, 4)
+min <- rep(0.1, 4)
+max <- rep(0.6, 4)
+toc <- 0.8
+lambda <- 0.25
+mu <- colMeans(R)
+target <- 0.0071
+
+group_mat <- cbind(c(1, 1, 0, 0), c(0, 0, 1, 1))
+grp_min <- c(0.05, 0.05)
+grp_max <- c(0.85, 0.85)
+
# Computes optimal weights for global minimum variance portfolio with
# constraints including turnover constraint
-turnover.gmv <- function(R, toc, weight.i, min, max){
-
- # number of assets in R
- p <- ncol(R)
-
- # Modify the returns matrix. This is done because there are 3 sets of
- # variables w.initial, w.buy, and w.sell
- returns <- cbind(R, R, R)
-
- V <- cov(returns)
- V <- make.positive.definite(V)
-
- # A matrix for initial weights
- A2 <- cbind(rep(1, p*3), rbind(diag(p), matrix(0, nrow=2*p, ncol=p)))
-
- # A matrix for lower box constraints
- Alo <- rbind(diag(p), diag(p), diag(p))
-
- # A matrix for upper box constraints
- Aup <- rbind(-diag(p), -diag(p), -diag(p))
-
- # vector to appply turnover constraint
- A3 <- c(rep(0, p), rep(-1, p), rep(1, p))
-
- # matrix for positive weight
- A4 <- rbind(matrix(0, nrow=p, ncol=p), diag(p), matrix(0, nrow=p, ncol=p))
-
- # matrix for negative weight
- A5 <- rbind(matrix(0, nrow=p*2, ncol=p), -diag(p))
-
- # Combine the temporary A matrices
- A.c <- cbind(A2, Alo, Aup, A3, A4, A5)
-
- # b vector holding the values of the constraints
- b <- c(1, weight.i, min, -max, -toc, rep(0, 2*p))
-
- # no linear term so set this equal to 0s
- d <- rep(0, p*3)
-
- sol <- solve.QP(Dmat=V, dvec=d, Amat=A.c, bvec=b, meq=6)
- wts <- sol$solution
- wts.final <- wts[(1:p)] + wts[(1+p):(2*p)] + wts[(2*p+1):(3*p)]
- wts.final
-}
-data(edhec)
-ret <- edhec[,1:5]
+# number of assets in R
+p <- ncol(R)
-# box constraints min and max
-min <- rep(0.1, 5)
-max <- rep(0.6, 5)
+# Modify the returns matrix. This is done because there are 3 sets of
+# variables w.initial, w.buy, and w.sell
+returns <- cbind(R, R, R)
-# turnover constraint
-toc <- 0.3
+V <- cov(returns)
+V <- make.positive.definite(V)
-# Initial weights vector
-weight.i <- rep(1/5,5)
+# A matrix for full investment, mean, and initial weights
+A2 <- cbind(rbind(diag(p), matrix(0, nrow=2*p, ncol=p)),
+ rep(mu, 3),
+ rep(1, p*3), rep(-1, p*3))
-opt.wts <- turnover.gmv(R=ret, toc=toc, weight.i=weight.i, min=min, max=max)
-opt.wts
+# A matrix for lower box constraints
+Alo <- rbind(diag(p), diag(p), diag(p))
+# A matrix for upper box constraints
+Aup <- rbind(-diag(p), -diag(p), -diag(p))
+
+# vector to appply turnover constraint
+A3 <- c(rep(0, p), rep(-1, p), rep(1, p))
+
+# matrix for positive weight
+A4 <- rbind(matrix(0, nrow=p, ncol=p), diag(p), matrix(0, nrow=p, ncol=p))
+
+# matrix for negative weight
+A5 <- rbind(matrix(0, nrow=p*2, ncol=p), -diag(p))
+
+# Combine the temporary A matrices
+A.c <- cbind(A2, Alo, Aup, A3, A4, A5)
+
+# Constraints matrix for group_mat
+A.c <- cbind(A.c, rbind(group_mat, group_mat, group_mat))
+A.c <- cbind(A.c, rbind(-group_mat, -group_mat, -group_mat))
+
+# b vector holding the values of the constraints
+b <- c(init.weights, target, 0.99, -1.01, min, -max, -toc, rep(0, 2*p))
+b <- c(b, grp_min, -grp_max)
+
+# no linear term so set this equal to 0s
+d <- rep(0, p*3)
+# d <- rep(-mu, 3)
+
+sol <- try(solve.QP(Dmat=make.positive.definite(2*lambda*V), dvec=d,
+ Amat=A.c, bvec=b, meq=4), silent=TRUE)
+if(inherits(sol, "try-error")) message("No solution found")
+wts <- sol$solution
+wts.final <- wts[(1:p)] + wts[(1+p):(2*p)] + wts[(2*p+1):(3*p)]
+wts.final
+sum(wts.final)
+wts.final %*% mu
+
# calculate turnover
-to <- sum(abs(diff(rbind(weight.i, opt.wts))))
+to <- sum(abs(diff(rbind(init.weights, wts.final))))
to
##### ROI Turnover constraints using ROI solver #####
-# Not working correctly.
-# Getting a solution now, but results are different than turnover.gmv
-# library(ROI)
-# library(ROI.plugin.quadprog)
-#
-#
-# # Use the first 5 funds in edhec for the returns data
-# ret <- edhec[, 1:5]
-# returns <- cbind(ret, ret, ret)
-#
-# V <- cov(returns)
-# V <- corpcor:::make.positive.definite(V)
-# mu <- apply(returns, 2, mean)
-# # number of assets
-# N <- ncol(returns)
-#
-# # Set the box constraints for the minimum and maximum asset weights
-# min <- rep(0.1, N/3)
-# max <- rep(0.6, N/3)
-#
-# # Set the bounds
-# bnds <- list(lower = list(ind = seq.int(1L, N/3), val = as.numeric(min)),
-# upper = list(ind = seq.int(1L, N/3), val = as.numeric(max)))
-# lambda <- 1
-# ROI_objective <- ROI:::Q_objective(Q=2*lambda*V, L=-mu*0)
-#
-# # Set up the Amat
-# # min_sum and max_sum of weights
-# A1 <- rbind(rep(1, N), rep(1, N))
-#
-# # initial weight matrix
-# A.iw <- cbind(diag(N/3), matrix(0, nrow=N/3, ncol=2*N/3))
-#
-# # turnover vector
-# A.t <- c(rep(0, N/3), rep(-1, N/3), rep(1, N/3))
-#
-# A.wpos <- t(cbind(rbind(matrix(0, ncol=N/3, nrow=N/3), diag(N/3), matrix(0, ncol=N/3, nrow=N/3)),
-# rbind(matrix(0, ncol=N/3, nrow=2*N/3), -diag(N/3))))
-#
-# Amat <- rbind(A1, A.iw, A.t, A.wpos)
-#
-# dir.vec <- c(">=","<=", rep("==", N/3), "<=", rep(">=", 2*N/3))
-# min_sum=1
-# max_sum=1
-# w.init <- rep(1/5, 5)
-# toc <- 0.3
-# rhs.vec <- c(min_sum, max_sum, w.init, toc, rep(0, 2*N/3))
-#
-# opt.prob <- ROI:::OP(objective=ROI_objective,
-# constraints=ROI:::L_constraint(L=Amat, dir=dir.vec, rhs=rhs.vec),
-# bounds=bnds)
-# roi.result <- ROI:::ROI_solve(x=opt.prob, solver="quadprog")
-# wts.tmp <- roi.result$solution
-# wts <- wts.tmp[1:(N/3)] + wts.tmp[(N/3+1):(2*N/3)] + wts.tmp[(2*N/3+1):N]
-# wts
\ No newline at end of file
+N <- ncol(R)
+
+tmpR <- cbind(R, R, R)
+V <- var(tmpR)
+V <- corpcor:::make.positive.definite(V)
+# lambda <- 0.25
+# mu <- colMeans(R)
+
+# Amat for full investment constraint
+Amat <- rbind(rep(1, N*3), rep(1, N*3))
+rhs <- c(1, 1)
+dir <- c("==", "==")
+# dir <- c(">=", "<=")
+
+# Amat for initial weights
+Amat <- rbind(Amat, cbind(diag(N), matrix(0, nrow=N, ncol=N*2)))
+rhs <- c(rhs, init.weights)
+dir <- c(dir, rep("==", N))
+
+# Amat for lower box constraints
+Amat <- rbind(Amat, cbind(diag(N), diag(N), diag(N)))
+rhs <- c(rhs, min)
+dir <- c(dir, rep(">=", N))
+
+# Amat for upper box constraints
+Amat <- rbind(Amat, cbind(-diag(N), -diag(N), -diag(N)))
+rhs <- c(rhs, -max)
+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, -toc)
+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))
+
+# set up the quadratic objective
+ROI_objective <- Q_objective(Q=make.positive.definite(2 * lambda * V), L=rep(0, N*3))
+
+opt.prob <- OP(objective=ROI_objective,
+ constraints=L_constraint(L=Amat, dir=dir, rhs=rhs))
+roi.result <- ROI_solve(x=opt.prob, solver="quadprog")
+# not sure why no solution with ROI
+print.default(roi.result)
+
+# check that the same constraints matrix and rhs vectors are used
+all.equal(t(Amat), A.c, check.attributes=FALSE)
+all.equal(rhs, b)
+
+# run solve.QP using Amat and rhs from ROI problem
+qp.result <- solve.QP(Dmat=make.positive.definite(2*lambda*V),
+ dvec=rep(0, N*3), Amat=t(Amat), bvec=rhs, meq=6)
+
+# results with solve.QP are working, but not with ROI
+all.equal(qp.result$solution, sol$solution)
+all.equal(roi.result$solution, sol$solution)
+
+roi.result$solution
+qp.result$solution
+sol$solution
+
+
+
More information about the Returnanalytics-commits
mailing list