[Returnanalytics-commits] r2221 - in pkg/PortfolioAnalytics: R sandbox
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Aug 4 03:36:03 CEST 2012
Author: hezkyvaron
Date: 2012-08-04 03:36:01 +0200 (Sat, 04 Aug 2012)
New Revision: 2221
Added:
pkg/PortfolioAnalytics/sandbox/sample_pso.R
Modified:
pkg/PortfolioAnalytics/R/optimize.portfolio.R
pkg/PortfolioAnalytics/sandbox/testing_ROI.R
Log:
- fixed some of the pso code in optimize.portfolio
- updated testing_ROI to reflect new integrated ROI
- wrote a similar test script for pso with comparison with DEoptim
- added a pso/DEoptim sample script to check speed and solution of both solvers, approach used is from Ardia et al., R Journal June 2011.
Modified: pkg/PortfolioAnalytics/R/optimize.portfolio.R
===================================================================
--- pkg/PortfolioAnalytics/R/optimize.portfolio.R 2012-08-03 10:47:37 UTC (rev 2220)
+++ pkg/PortfolioAnalytics/R/optimize.portfolio.R 2012-08-04 01:36:01 UTC (rev 2221)
@@ -265,41 +265,6 @@
-## The first draft of the integrated ROI:
-# if (optimize_method == "ROI") {
-# bnds <- list(lower = list(ind = seq.int(1L, N), val = constraints$min),
-# upper = list(ind = seq.int(1L, N), val = constraints$max))
-# objectives <- do.call(cbind, sapply(constraints$objectives,"[", "name"))
-# moments <- list()
-# for (i in 1:length(objectives)) moments[[i]] <- eval(as.symbol(objectives[i]))(R)
-# names(moments) <- objectives
-# plugin <- ifelse(any(objectives == "var"), "quadprog", "glpk")
-# target.return <- do.call(cbind, sapply(constraints$objectives,"[", "target"))
-# lambda <- do.call(cbind, sapply(constraints$objectives, "[", "risk_aversion"))
-# if (is.null(lambda)) lambda <- 1
-# if (plugin == "quadprog") ROI_objective <- ROI:::Q_objective(Q = 2 * lambda * moments$var, L = -moments$mean)
-# if (plugin == "glpk") ROI_objective <- ROI:::L_objective(L = -moments$mean)
-# Amat <- rbind(rep(1, N), rep(1, N))
-# dir.vec <- c(">=", "<=")
-# rhs.vec <- c(constraints$min_sum, constraints$max_sum)
-# if (!is.null(target.return)) {
-# Amat <- rbind(Amat, moments$mean)
-# dir.vec <- cbind(dir.vec, "==")
-# rhs.vec <- cbind(rhs.vec, target.return)
-# }
-# q.prob <- ROI:::OP(objective = ROI_objective,
-# constraints = L_constraint(L = Amat, dir = dir.vec, rhs = rhs.vec),
-# bounds = bnds)
-# roi.result <- ROI:::ROI_solve(x = q.prob, solver = plugin)
-# weights <- roi.result$solution
-# names(weights) <- colnames(R)
-# out$weights <- weights
-# out$objective_measures <- roi.result$objval
-# out$call <- call
-# }
-
-
-
if(optimize_method == "ROI"){
# This takes in a regular constraint object and extracts the desired business objectives
# and converts them to matrix form to be inputed into a closed form solver
@@ -309,18 +274,15 @@
# retrive the objectives to minimize, these should either be "var" and/or "mean"
# we can eight miniminze variance or maximize quiadratic utility (we will be minimizing the neg. quad. utility)
moments <- list(mean=rep(0, N), var=NULL)
- target <- NULL
- lambda <- NULL
alpha <- 0.05
for(objective in constraints$objectives){
if(objective$enabled){
- if(objective$name != "mean" || objective$name != "var" || objective$name != "CVaR")
- stop("ROI only solves mean, var, or sample CVaR type business objectives,
- choose a different optimize_method.")
+ if(!any(c(objective$name == "mean", objective$name == "var", objective$name == "CVaR")))
+ stop("ROI only solves mean, var, or sample CVaR type business objectives, choose a different optimize_method.")
moments[[objective$name]] <- eval(as.symbol(objective$name))(R)
target <- ifelse(!is.null(objective$target),objective$target, NA)
alpha <- ifelse(!is.null(objective$alpha), objective$alpha, alpha)
- lambda <- ifelse(!is.null(objective$risk_aversion),objective$risk_aversion,1)
+ lambda <- ifelse(!is.null(objective$risk_aversion), objective$risk_aversion, 1)
}
}
plugin <- ifelse(any(names(moments)=="var"), "quadprog", "glpk")
@@ -356,9 +318,8 @@
## case if method=pso---particle swarm
if(optimize_method=="pso"){
stopifnot("package:pso" %in% search() || require("pso",quietly = TRUE) )
- # DEoptim does 200 generations by default, so lets set the size of each generation to search_size/200)
if(hasArg(maxit)) maxit=match.call(expand.dots=TRUE)$maxit else maxit=N*50
- PSOcformals <- formals(psoptim.control)
+ PSOcformals <- list(trace=FALSE, fnscale=1, maxit=1000, maxf=Inf, abstol=-Inf, reltol=0)
PSOcargs <- names(PSOcformals)
if( is.list(dotargs) ){
pm <- pmatch(names(dotargs), PSOcargs, nomatch = 0L)
Added: pkg/PortfolioAnalytics/sandbox/sample_pso.R
===================================================================
--- pkg/PortfolioAnalytics/sandbox/sample_pso.R (rev 0)
+++ pkg/PortfolioAnalytics/sandbox/sample_pso.R 2012-08-04 01:36:01 UTC (rev 2221)
@@ -0,0 +1,37 @@
+# # # # # # # # # #
+# Try DEoptim
+library(xts)
+library(quantmod)
+library(PerformanceAnalytics)
+library(DEoptim)
+
+data(edhec)
+R <- edhec
+N <- ncol(edhec)
+mu <- colMeans(R)
+sigma <- cov(R)
+
+obj <- function(w){
+ if (sum(w)==0) {w <- w + 1e-2}
+ w <- w/sum(w)
+ CVaR <- ES(weights= w,
+ method="gaussian",
+ portfolio_method="component",
+ mu=mu,
+ sigma=sigma)
+ tmp1 <- CVaR$ES
+ tmp2 <- max(CVaR$pct_contrib_ES - 0.225, 0)
+ out <- tmp1 + 1e3 * tmp2
+}
+
+set.seed(1234)
+out <- DEoptim(fn=obj,
+ lower = rep(0, N),
+ upper = rep(1, N))
+
+out$optim$bestval
+wts.deoptim <- out$optim$bestmem / sum(out$optim$bestmem)
+
+test <- psoptim(rep(NA,N), obj, lower=0, upper=5, control=list(abstol=1e-8))
+test$value
+wts.pso <- test$par/sum(test$par)
Modified: pkg/PortfolioAnalytics/sandbox/testing_ROI.R
===================================================================
--- pkg/PortfolioAnalytics/sandbox/testing_ROI.R 2012-08-03 10:47:37 UTC (rev 2220)
+++ pkg/PortfolioAnalytics/sandbox/testing_ROI.R 2012-08-04 01:36:01 UTC (rev 2221)
@@ -1,5 +1,5 @@
-# # # # # # # # # # # # # #
-# OPTIMIZATION TESTING
+# # # # # # # # # # # # # # # # #
+# OPTIMIZATION TESTING: ROI
#
library(xts)
@@ -12,191 +12,63 @@
library(Ecdat)
library(PortfolioAnalytics)
+
+# General Parameters for sample code
data(edhec)
-cov.mat <- var(edhec)
-mu.vec <- apply(edhec, 2, mean)
-n.assets <- ncol(edhec)
+funds <- names(edhec)
+mu.port <- mean(colMeans(edhec))
+gen.constr <- constraint(assets = colnames(edhec), min=-Inf, max =Inf, min_sum=1, max_sum=1, risk_aversion=1)
+gen.constr <- add.objective(constraints=no.box.constr, type="return", name="mean", enabled=FALSE, multiplier=0, target=mu.port)
+gen.constr <- add.objective(constraints=no.box.constr, type="risk", name="var", enabled=FALSE, multiplier=0, risk_aversion=10)
+gen.constr <- add.objective(constraints=no.box.constr, type="risk", name="CVaR", enabled=FALSE, multiplier=0)
-# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
-# Sample portfolio optimization problems....
-
# =====================
-# Arbitrage
-#
-# to check the values of this example, I did the same problem with Rnupot.
-# unfortunately, ROI (or Rglpk, more specifically) cannot solve this problem well.
-# when using a RHS term in the constraint of Ax >= 0, the solution is trivial,
-# with RHS = rep(0.001, n.assets), the solution is non-trivial, but wrong.
+# Max return under box constraints, fully invested
#
-# data(edhec)
-# set.seed(123)
-# n.assets <- ncol(edhec)
-# S <- 1 + apply(edhec, 2, sample, n.assets)
-# solution <- solveQP(objL=rep(1, ncol(S)), S,
-# cLO=rep(0, nrow(S)), cUP=rep(Inf, nrow(S)),
-# type=minimize)
-# wts <- solution$variables$x$current
-# > wts
-# 1 2 3 4 5 6 7
-# -0.03687089 0.34978850 -0.02893078 0.15988671 0.29418430 0.02867836 -0.29697380
-# 8 9 10 11 12 13
-# -0.45799863 -0.30662611 -0.11010508 -0.74430483 -0.02668694 0.18126218
+max.port <- gen.constr
+max.port$min <- 0.01
+max.port$max <- 0.30
+max.port$objectives[[1]]$enabled <- TRUE
+max.port$objectives[[1]]$target <- NULL
+max.solution <- optimize.portfolio(edhec, max.port, "ROI")
-set.seed(123)
-S.mat <- 1 + apply(edhec, 2, sample, size=n.assets)
-bnds <- list(lower = list(ind = seq.int(1L, n.assets), val = rep(-Inf,n.assets)),
- upper = list(ind = seq.int(1L, n.assets), val = rep(Inf,n.assets)))
-arb.prob <- OP(objective = L_objective(L=rep(1, n.assets)),
- constraints = L_constraint(L=S.mat,
- dir=rep(">=", n.assets),
- rhs=rep(0.001, n.assets)),
- bounds=bnds)
-arb.constr <- constraint_ROI(op.problem=arb.prob, "glpk")
-arb.test <- optimize.portfolio(edhec, arb.constr, "ROI")
-# > test$solution with ROI, same is found via optimize.portfolio(..., method="ROI")
-# [1] -0.37584904 0.00000000 0.00000000 0.18434868 0.00000000 0.00000000 0.00000000 0.32167429 0.14030506
-# [10] 0.00000000 -0.41200277 -0.02153654 0.25916374
-#
-# Issues that I have found:
-# 1) arb.test$constraints returns:
-#
-# > test$constraints
-# An object containing 6 constraints.
-# Some constraints are of type nonlinear.
-#
-# which is meaningless,
-#
-# 2) when not passing in assets into arb.prob, this is returned:
-#
-# assuming equal weighted seed portfolio
-# assuming equal weighted seed portfolio
-#
-# only should say it once.
-
-
# =====================
-# Mean-variance: Maximize quadratic utility
+# Mean-variance: Fully invested, Global Minimum Variance Portfolio
#
-funds <- names(edhec)
-mu.port <- 0.002
-Amat <- rbind(rep(1,n.assets), mu.vec)
-lambda <- 1
-mean.var.prob <- OP(objective=Q_objective(Q=2*lambda*cov.mat, L=-mu.vec),
- constraints=L_constraint(L=Amat,
- dir=c("==","=="),
- rhs=c(1,mu.port)),
- bounds=bnds)
-mean.var.constr <- constraint_ROI(assets=funds, op.problem=mean.var.prob, solver="quadprog")
-wts <- ROI_solve(x=mean.var.prob, solver="quadprog")$solution
-mean.var.solution <- optimize.portfolio(edhec, mean.var.constr, "ROI_old")
-# using integrated ROI
+gmv.port <- gen.constr
+gmv.port$objectives[[2]]$enabled <- TRUE
+gmv.port$objectives[[2]]$risk_aversion <- 1
+gmv.solution <- optimize.portfolio(edhec, gmv.port, "ROI")
-mean.var <- constraint(assets = colnames(edhec), min=-Inf, max =Inf, min_sum=1, max_sum=1, risk_aversion=1)
-mean.var <- add.objective(constraints=mean.var, type="return", name="mean", enabled=TRUE, multiplier=0, target=mu.port)
-mean.var <- add.objective(constraints=mean.var, type="risk", name="var", enabled=TRUE, multiplier=0)
-solution <- optimize.portfolio(edhec, mean.var, "ROI")
-paste(names(edhec),solution$weights)
-# results for this are:
+# ========================
+# Mean-variance: Maximize quadratic utility, fully invested, target portfolio return
#
-# > mean.var.solution$weights
-# Convertible Arbitrage CTA Global Distressed Securities Emerging Markets
-# -0.38704286 0.08062104 -0.35410876 -0.06088908
-# Equity Market Neutral Event Driven Fixed Income Arbitrage Global Macro
-# 0.49829093 -0.49045547 0.74727967 -0.49173927
-# Long/Short Equity Merger Arbitrage Relative Value Short Selling
-# -0.44361180 0.52390274 0.37357962 -0.04284137
-# Funds of Funds
-# 1.04701463
-# > wts
-# [1] -0.38704286 0.08062104 -0.35410876 -0.06088908 0.49829093 -0.49045547 0.74727967
-# [8] -0.49173927 -0.44361180 0.52390274 0.37357962 -0.04284137 1.04701463
-#
-# however, the constraints show up as:
-# $constraints
-# An object containing 6 constraints.
-# Some constraints are of type nonlinear.
+target.port <- gen.constr
+target.port$objectives[[1]]$enabled <- TRUE
+target.port$objectives[[2]]$enabled <- TRUE
+target.solution <- optimize.portfolio(edhec, target.port, "ROI")
-
-# ========================================================
-# Mean-variance: Maximize quadratic utility --- dollar neutral
+# ========================
+# Mean-variance: Maximize quadratic utility, dollar-neutral, target portfolio return
#
-funds <- names(edhec)
-mu.port <- 0.002
-Amat <- cbind(rep(1,n.assets), mu.vec)
-dollar.neutral.prob <- OP(objective=Q_objective(Q=2*cov.mat, L=-mu.vec),
- constraints=L_constraint(L=t(Amat),
- dir=c("==","=="),
- rhs=c(0,mu.port)),
- bounds=bnds)
-dollar.neutral.constr <- constraint_ROI(assets=funds, op.problem=dollar.neutral.prob, solver="quadprog")
-wts <- ROI_solve(x=dollar.neutral.prob, solver="quadprog")$solution
-dollar.neutral.solution <- optimize.portfolio(edhec, dollar.neutral.constr, "ROI")
-paste(funds, dollar.neutral.solution$weights)
+dollar.neu.port <- gen.constr
+dollar.neu.port$min_sum <- 0
+dollar.neu.port$max_sum <- 0
+dollar.neu.port$objectives[[1]]$enabled <- TRUE
+dollar.neu.port$objectives[[2]]$enabled <- TRUE
+dollar.neu.solution <- optimize.portfolio(edhec, dollar.neu.port, "ROI")
-# using integrated ROI
-mean.var <- constraint(assets = colnames(edhec), min = -Inf, max = Inf, min_sum=0, max_sum=0, risk_aversion=1)
-mean.var <- add.objective(constraints=mean.var, type="return", name="mean", enabled=TRUE, multiplier=0, target=mu.port)
-mean.var <- add.objective(constraints=mean.var, type="risk", name="var", enabled=TRUE, multiplier=0)
-solution <- optimize.portfolio(edhec, mean.var, "ROI_new")
-paste(names(edhec),solution$weights)
-
-# =====================
-# Maximize return given box constraints
+# ========================
+# Minimize CVaR with target return
#
-# A set of box constraints used to initialize portfolios
-init.constr <- constraint(assets = colnames(edhec), min = .05, max = .3, min_sum=1, max_sum=1)
-init.constr <- add.objective(constraints=init.constr, type="return", name="mean", enabled=TRUE, multiplier=0)
-test <- optimize.portfolio(edhec, init.constr, "ROI_new")
+cvar.port <- gen.constr
+cvar.port$objectives[[1]]$enabled <- TRUE
+cvar.port$objectives[[3]]$enabled <- TRUE
+cvar.solution <- optimize.portfolio(edhec, cvar.port, "ROI")
-
-
-
-# =====================
-# Mean-variance: Maximize return, with constraint on variance
-# setting varaince to be the average variance between the edhec funds.
-#
-avg.var <- mean(apply(edhec, 2, var))
-max.mean.prob <- OP(objective=L_objective(mu.vec),
- constraints=Q_constraint(Q=list(2*cov.mat, matrix(0, nrow=n.assets, ncol=n.assets)),
- L=rbind(rep(0, n.assets),rep(1,n.assets)),
- dir=c("==","=="),
- rhs=c(avg.var,1)),
- bounds=bnds,
- maximum=TRUE)
-max.mean.constr <- constraint_ROI(assets=funds, op.problem=max.mean.prob, solver="glpk")
-wts <- ROI_solve(x=max.mean.prob, solver="glpk")$solution
-max.mean.solution <- optimize.portfolio(edhec, max.mean.constr, "ROI")
-#
-# As expected, attempting to set up this problem leads to failure.
-# need to implement a second order conic solver.
-# this solver is in the package CLSOCP
-
-
-# ===================
-# This secion works, and I am puting it aside for now
-#
-# # Comparing resutls wtih Guy's slides of PortfolioOptimization
-# # sllide number 24/70, mean-variance optimization
-# # subject to fully-invested and expected portfolio return constraints
-# data(CRSPday)
-# R <- 100*CRSPday[,4:6]
-# mean_vect <- apply(R,2,mean)
-# cov_mat <- var(R)
-# Amat <- rbind(rep(1,3),mean_vect)
-# mu.port <- 0.1
-# bnds <- list(lower = list(ind = seq.int(1L, as.integer(3)), val = rep(-Inf,3)),
-# upper = list(ind = seq.int(1L, as.integer(3)), val = rep(Inf,3)))
-# q.prob <- OP(objective=Q_objective(Q=2*cov_mat, L=rep(0,3)),
-# constraints=L_constraint(L=Amat,
-# dir=c("==","=="),
-# rhs=c(1, mu.port)),
-# bounds=bnds)
-# wts <- ROI_solve(x=q.prob, solver="quadprog")$solution
-
-
More information about the Returnanalytics-commits
mailing list