R/methods-PMSOLVER.R

In pmoments: Partial Moments

#################################################################################
##
##   R package pmoments by Alexios Ghalanos Copyright (C) 2008
##   This file is part of the R package pmoments.
##
##   The R package pmoments is free software: you can redistribute it and/or modify
##   it under the terms of the GNU General Public License as published by
##   the Free Software Foundation, either version 3 of the License, or
##   (at your option) any later version.
##
##   The R package pmoments is distributed in the hope that it will be useful,
##   but WITHOUT ANY WARRANTY; without even the implied warranty of
##   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##   GNU General Public License for more details.
##
#################################################################################

pmSolver<-function(pcmExpectation, userConstraints, npoints=200, type=c("I","II","III","IV"))
{
	UseMethod("pmSolver")
}

.pmSolver<-function(pcmExpectation, userConstraints, npoints=200, type=c("I","II","III","IV"))
{
	if(length(type)!=1) type = "I"
	sol<-switch(type,
			I = .pmSolverI(pcmExpectation, userConstraints),
			II = .pmSolverII(pcmExpectation, userConstraints),
			III = .pmSolverIII(pcmExpectation, userConstraints),
			IV = .pmSolverIV(pcmExpectation, userConstraints, npoints))
	return(sol)
	
}

setMethod("pmSolver", signature(pcmExpectation="PCME",userConstraints="userConstraints"), .pmSolver)

# Solver I : Minimum Risk Portfolio
.pmSolverI<-function(pcmExpectation, userConstraints)
{
	if(dim(pcmExpectation@pm)[2]!=length(userConstraints@forecasts)) stop("length of forecasts not equal to partial co-moment width")
	n<-length(userConstraints@forecasts)
	# we use the quadratic solver of R which has an L00 representation with direction = -1 (solverConstraints methodology).
	scon<-solverConstraints(userConstraints, type = "L00", direction = -1)
	lambda<-userConstraints@riskAversion
	m = pcmExpectation@moment
	if(m==0) mm=1 else mm=m
	if(m==0) mm=1 else mm=m
	slpm<-pcmExpectation@pm
	riskFree=userConstraints@riskFree
	# Last 2 lines of bvec/Amat are the forecasts
	bvec<-scon@bineq
	Amat<-scon@Aineq
	nn=dim(Amat)[1]
	bvec = bvec[-c(nn-1,nn)]
	Amat = Amat[-c(nn-1,nn),]
	dvec<-rep(0,n)
	# implement try catch and create solverMessage
	xtmp<-try(sol<-solve.QP(Dmat=(lambda*slpm), dvec=dvec, Amat=t(Amat), bvec=bvec),silent = TRUE)
	if(inherits(xtmp, "try-error")){
		smessage="non-convergence"
		reward = NA
		w=matrix(rep(NA,n),ncol=n)
		risk=NA
	}
	else{
		smessage="convergence"
		w = matrix(sol$solution,ncol=n)
		risk = (sol$solution%*%slpm%*%sol$solution)^(1/mm)
		reward = userConstraints@forecasts%*%as.numeric(sol$solution)
	}
	if(!is.null(names(userConstraints@forecasts))) assetnames = names(userConstraints@forecasts) else assetnames=paste(rep("Asset ",n),1:n,sep="")
	colnames(w)<-assetnames
	ans=new("PMSOLVER",
			weights = w,
			riskMeasure = risk,
			rewardMeasure = reward,
			riskAversion = lambda,
			riskFree = riskFree,
			tangent = NA,
			moment=m,
			solverMessage = smessage,
			solver = "quadratic",
			solverType = "Partial Moments (I)",
			assetnames = colnames(pcmExpectation@pm))
	return(ans)
}
# Method 2: Minimize risk for a given level of return
.pmSolverII<-function(pcmExpectation, userConstraints)
{
	# check that constraints and pcmExpectation are same dimensions in width
	if(dim(pcmExpectation@pm)[2]!=length(userConstraints@forecasts)) stop("length of forecasts not equal to partial co-moment width")
	n<-length(userConstraints@forecasts)
	# we use the quadratic solver of R which has an L00 representation with direction = -1 (solverConstraints methodology).
	scon<-solverConstraints(userConstraints, type = "L00", direction = -1)
	lambda<-userConstraints@riskAversion
	m = pcmExpectation@moment
	if(m==0) mm=1 else mm=m
	slpm<-pcmExpectation@pm
	riskFree=userConstraints@riskFree
	bvec<-scon@bineq
	Amat<-scon@Aineq
	dvec<-rep(0,n)
	# implement try catch and create solverMessage
	xtmp<-try(sol<-solve.QP(Dmat=(lambda*slpm), dvec=dvec, Amat=t(Amat), bvec=bvec),silent = TRUE)
	if(inherits(xtmp, "try-error")){
		smessage="non-convergence"
		reward = NA
		w=matrix(rep(NA,n),ncol=n)
		risk=NA
	}
	else{
		smessage="convergence"
		w = matrix(sol$solution,ncol=n)
		risk = (sol$solution%*%slpm%*%sol$solution)^(1/mm)
		reward = userConstraints@forecasts%*%as.numeric(sol$solution)
	}
	if(!is.null(names(userConstraints@forecasts))) assetnames = names(userConstraints@forecasts) else assetnames=paste(rep("Asset ",n),1:n,sep="")
	colnames(w)<-assetnames
	ans=new("PMSOLVER",
			weights = w,
			riskMeasure = risk,
			rewardMeasure = reward,
			riskAversion=lambda,
			riskFree = riskFree,
			tangent = NA,
			moment=m,
			solverMessage = smessage,
			solver = "quadratic",
			solverType = "Partial Moments (II)",
			assetnames = colnames(pcmExpectation@pm))
	return(ans)
}

# Method 3: Maximize the risk adjusted return for a given risk aversion coefficient
.pmSolverIII<-function(pcmExpectation, userConstraints)
{
	# check that constraints and pcmExpectation are same dimensions in width
	if(dim(pcmExpectation@pm)[2]!=length(userConstraints@forecasts)) stop("length of forecasts not equal to partial co-moment width")
	n<-length(userConstraints@forecasts)
	lambda<-userConstraints@riskAversion
	# we use the quadratic solver of R which has an L00 representation with direction = -1 (solverConstraints methodology).
	scon<-solverConstraints(userConstraints, type = "L00", direction = -1)
	m = pcmExpectation@moment
	if(m==0) mm=1 else mm=m
	slpm<-pcmExpectation@pm
	riskFree=userConstraints@riskFree
	# we need to remove forecast/return constraint from this method
	# this is easy as it is always the last line(s) in the matrix
	bvec<-scon@bineq[-c(length(scon@bineq),length(scon@bineq)-1)]
	Amat<-scon@Aineq[-c(dim(scon@Aineq)[1],dim(scon@Aineq)[1]-1),]
	dvec<-userConstraints@forecasts
	# implement try catch and create solverMessage
	# 2x lambda since already 0.5x lambda (standardization)
	xtmp<-try(sol<-solve.QP(Dmat=2*lambda*slpm, dvec=dvec, Amat=t(Amat), bvec=bvec),silent = TRUE)
	if(inherits(xtmp, "try-error")){
		smessage="non-convergence"
		w=matrix(rep(NA,n),ncol=n)
		risk=NA
		reward=NA
	}
	else{
		w = matrix(sol$solution,ncol=n)
		risk<-(sol$solution%*%slpm%*%sol$solution)^(1/mm)
		reward<-userConstraints@forecasts%*%sol$solution
		smessage="convergence"
	}
	if(!is.null(names(userConstraints@forecasts))) assetnames = names(userConstraints@forecasts) else assetnames=paste(rep("Asset ",n),1:n,sep="")
	colnames(w)<-assetnames
	ans<-new("PMSOLVER",
			weights = w,
			riskMeasure = risk,
			rewardMeasure = reward,
			riskAversion=lambda,
			riskFree = riskFree,
			tangent = NA,
			moment=m,
			solverMessage = smessage,
			solver = "quadratic",
			solverType = "Partial Moments (III)",
			assetnames = colnames(pcmExpectation@pm))
	return(ans)
}

# Method 4: PM Frontier
.pmSolverIV<-function(pcmExpectation, userConstraints, npoints=200)
{
	if(dim(pcmExpectation@pm)[2]!=length(userConstraints@forecasts)) stop("length of forecasts not equal to partial co-moment width")
	n<-length(userConstraints@forecasts)
	# lambd: 1 is the cuttof for optimality, 2 makes it look nice by bending backwards into non optimal territory
	lambda<-seq(0,2,length.out=npoints)
	# we use the quadratic solver of R which has an L00 representation with direction = -1 (solverConstraints methodology).
	scon<-solverConstraints(userConstraints, type = "L00", direction = -1)
	m = pcmExpectation@moment[1]
	slpm<-pcmExpectation@pm
	riskFree=userConstraints@riskFree
	if(m==0) mm=1 else mm=m
	# we need to remove forecast/return constraint from this method
	# this is easy as it is always the last line(s) in the matrix
	bvec<-scon@bineq[-c(length(scon@bineq),length(scon@bineq)-1)]
	Amat<-scon@Aineq[-c(dim(scon@Aineq)[1],dim(scon@Aineq)[1]-1),]
	dvec<-userConstraints@forecasts
	# implement try catch and create solverMessage
	w<-matrix(NA,ncol=n,nrow=npoints)
	risk<-vector(mode="numeric",length=npoints)
	reward<-vector(mode="numeric",length=npoints)
	smessage<-vector(mode="character",length=npoints)
	for(i in 1:npoints){
		xtmp<-try(sol<-solve.QP(Dmat=(2*lambda[i])*slpm, dvec=(1-lambda[i])*dvec, Amat=t(Amat), bvec=bvec), silent = TRUE)
		if(inherits(xtmp, "try-error")){
			smessage[i]="non-convergence"
			w[i,]=rep(NA,n)
			risk[i]=NA
			reward[i]=NA
		}
		else{
			w[i,] = sol$solution
			risk[i]<-(sol$solution%*%slpm%*%sol$solution)^(1/mm)
			reward[i]<-as.numeric(userConstraints@forecasts)%*%sol$solution
			smessage[i]="convergence"
		}
	}
	
	valid<-which(!is.na(risk))
	tangent = which(max((reward[valid]-riskFree)/risk[valid]) == (reward[valid]-riskFree)/risk[valid])
	if(!is.null(names(userConstraints@forecasts))) assetnames = names(userConstraints@forecasts) else assetnames=paste(rep("Asset ",n),1:n,sep="")
	colnames(w)<-assetnames
	ans<-new("PMSOLVER",
			weights = w[valid,],
			riskMeasure = risk[valid],
			rewardMeasure = reward[valid],
			riskAversion=lambda,
			riskFree = riskFree,
			tangent = tangent,
			moment=m,
			solverMessage = smessage,
			solver = "quadratic",
			solverType = "Partial Moments (IV)",
			assetnames = colnames(pcmExpectation@pm))
	return(ans)
}

.getweights<-function(x, row.names = NULL, optional = FALSE, ....)
{
	as.data.frame(x@weights)
}
setMethod("as.data.frame",signature(x="PMSOLVER"),.getweights)