parma
Portfolio Allocation and Risk Management Applications

Package index
Search the parma package
Vignettes
Functions
314
Source code
17
Man pages
15
Home
/
CRAN
/
parma
/
R/p-MILP.R

R/p-MILP.R
In parma: Portfolio Allocation and Risk Management Applications

Defines functions milpfracsolver milpminsolver milpport milpcdar.minrisk milpcvar.minrisk milpminmax.minrisk milpmad.minrisk milplpm1.minrisk milpmodrachev 

#################################################################################
##
##   R package parma
##   Alexios Galanos Copyright (C) 2012-2013 (<=Aug)
##   Alexios Galanos and Bernhard Pfaff Copyright (C) 2013- (>Aug)
##   This file is part of the R package parma.
##
##   The R package parma 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 parma 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.
##
#################################################################################


#--------------------------------------------------------------------------------
# Rachev Ratio: Upper to Lower CVaR
# Research Reference (Rachev Ratio) : Biglova, A., Ortobelli, S., Rachev, S., Stoyanov (2004)
# Research Reference (Mixed Integer LP representation) : Konno, Tanaka, Yamamoto (2008)
#--------------------------------------------------------------------------------
# This expresses the fractional function as a MILP problem by applying the Charnes-Cooper Transformation
# Alexios Galanos 2010

# Completely unrealistic to expect to solve for scenario size>200 and alpha>0.01
# (at least with GLPK...cplex has more power but still very difficult).

milpmodrachev = function(optvars)
{
	Data = optvars$Data
	n = NROW(Data)
	m = NCOL(Data)
	colnames(Data) = NULL
	rownames(Data) = NULL
	probability = optvars$probability
	if(is.null(probability)) probability = rep(1/n, n)
	Data = scale(Data, scale = FALSE)
	
	upa = max(1, as.integer(optvars$upper.alpha * n))
	lpa = max(1, as.integer(optvars$lower.alpha * n))
	
	# We need specialized inequalities and parameter bounds inequalities multiplied by y0
	
	# S, alpha, theta, u_t, v_t, phi_t, psi_t, z_t, slack1_t, slack2_t, constant, y0
	
	# Constraint1: theta + sum(psi_t)/(1 - upper.alpha)*n = 1
	Amat = cbind(
			matrix(0, ncol = m, nrow = 1), 
			matrix(0, ncol = 1, nrow = 1), 
			matrix(1, ncol = 1, nrow = 1),
			matrix(0, ncol = n, nrow = 1),
			matrix(0, ncol = n, nrow = 1),
			matrix(0, ncol = n, nrow = 1),
			matrix(probability/upa, ncol = n, nrow = 1),
			matrix(0, ncol = n, nrow = 1),
			matrix(0, ncol = n, nrow = 1),
			matrix(0, ncol = n, nrow = 1),
			matrix(0, ncol = 1, nrow = 1),
			matrix(0, ncol = 1, nrow = 1)
	)
	
	# Constraint2: u_t - v_t + sum(w * S_t) + a = 0
	Amat = rbind( Amat, cbind( 
					Data,
					matrix(1, ncol = 1, nrow = n),
					matrix(0, ncol = 1, nrow = n),
					diag(n),
					-1 * diag(n),
					matrix(0, ncol = n, nrow = n),
					matrix(0, ncol = n, nrow = n),
					matrix(0, ncol = n, nrow = n),
					matrix(0, ncol = n, nrow = n),
					matrix(0, ncol = n, nrow = n),
					matrix(0, ncol = 1, nrow = n),
					matrix(0, ncol = 1, nrow = n)) )
	
	# Constraint3: phi_t - psi_t - sum(w * S_t) + theta = 0
	Amat = rbind( Amat, cbind( 
					-1 * Data,
					matrix(0, ncol = 1, nrow = n),
					matrix(1, ncol = 1, nrow = n),
					matrix(0, ncol = n, nrow = n),
					matrix(0, ncol = n, nrow = n),
					diag(n),
					-1 * diag(n),
					matrix(0, ncol = n, nrow = n),
					matrix(0, ncol = n, nrow = n),
					matrix(0, ncol = n, nrow = n),
					matrix(0, ncol = 1, nrow = n),
					matrix(0, ncol = 1, nrow = n)) )
	
	# Constraint4: phi_t - 1000*z_t + slack1_t = 0
	Amat = rbind( Amat, cbind( 
					matrix(0, ncol = m, nrow = n),
					matrix(0, ncol = 1, nrow = n),
					matrix(0, ncol = 1, nrow = n),
					matrix(0, ncol = n, nrow = n),
					matrix(0, ncol = n, nrow = n),
					diag(n),
					matrix(0, ncol = n, nrow = n),
					-1000 * diag(n),
					diag(n),
					matrix(0, ncol = n, nrow = n),
					matrix(0, ncol = 1, nrow = n),
					matrix(0, ncol = 1, nrow = n)) )
	
	# Constraint5: psi_t - 1000 + 1000*z_t + slack3_t = 0
	Amat = rbind( Amat, cbind( 
					matrix(0, ncol = m, nrow = n),
					matrix(0, ncol = 1, nrow = n),
					matrix(0, ncol = 1, nrow = n),
					matrix(0, ncol = n, nrow = n),
					matrix(0, ncol = n, nrow = n),
					matrix(0, ncol = n, nrow = n),
					diag(n),
					1000 * diag(n),
					matrix(0, ncol = n, nrow = n),
					diag(n),
					matrix(-1000, ncol = 1, nrow = n),
					matrix(0, ncol = 1, nrow = n) ) )
	
	# Constraint6: sum(z_t) = (1 - b) * n
	Amat = rbind( Amat, cbind( 
					matrix(0, ncol = m, nrow = 1),
					matrix(0, ncol = 1, nrow = 1),
					matrix(0, ncol = 1, nrow = 1),
					matrix(0, ncol = n, nrow = 1),
					matrix(0, ncol = n, nrow = 1),
					matrix(0, ncol = n, nrow = 1),
					matrix(0, ncol = n, nrow = 1),
					matrix(1, ncol = n, nrow = 1),
					matrix(0, ncol = n, nrow = 1),
					matrix(0, ncol = n, nrow = 1),
					matrix(0, ncol = 1, nrow = 1),
					matrix(0, ncol = 1, nrow = 1) ) )
	
	# Constraint7: sum(x_t) = budget*y0
	Amat = rbind( Amat, cbind(
					matrix(1, ncol = m, nrow = 1),
					matrix(0, ncol = 1, nrow = 1),
					matrix(0, ncol = 1, nrow = 1),
					matrix(0, ncol = n, nrow = 1),
					matrix(0, ncol = n, nrow = 1),
					matrix(0, ncol = n, nrow = 1),
					matrix(0, ncol = n, nrow = 1),
					matrix(0, ncol = n, nrow = 1),
					matrix(0, ncol = n, nrow = 1),
					matrix(0, ncol = n, nrow = 1),
					matrix(0, ncol = 1, nrow = 1),
					matrix(-optvars$budget, ncol = 1, nrow = 1) ) )
	
	# Constraint8: x<=UB.y0 i.e. x - UB.y0 <=0
	Amat = rbind( Amat, cbind( 
					diag(m),
					matrix(0, ncol = 1, nrow = m),
					matrix(0, ncol = 1, nrow = m),
					matrix(0, ncol = n, nrow = m),
					matrix(0, ncol = n, nrow = m),
					matrix(0, ncol = n, nrow = m),
					matrix(0, ncol = n, nrow = m),
					matrix(0, ncol = n, nrow = m),
					matrix(0, ncol = n, nrow = m),
					matrix(0, ncol = n, nrow = m),
					matrix(0, ncol = 1, nrow = m),
					matrix(-optvars$UB, ncol = 1, nrow = m) ) )
	
	# Constraint9: x - LB.y0 >=0
	Amat = rbind( Amat, cbind(
					diag(m),
					matrix(0, ncol = 1, nrow = m),
					matrix(0, ncol = 1, nrow = m),
					matrix(0, ncol = n, nrow = m),
					matrix(0, ncol = n, nrow = m),
					matrix(0, ncol = n, nrow = m),
					matrix(0, ncol = n, nrow = m),
					matrix(0, ncol = n, nrow = m),
					matrix(0, ncol = n, nrow = m),
					matrix(0, ncol = n, nrow = m),
					matrix(0, ncol = 1, nrow = m),
					matrix(-optvars$LB, ncol = 1, nrow = m) ) )
	
	# Constraint10: any user equalities
	if(!is.null(optvars$eq.mat)){
		eqm = dim(optvars$eq.mat)[1]
		Amat = rbind( Amat, cbind(
						optvars$eq.mat,
						matrix(0, ncol = 1, nrow = eqm),
						matrix(0, ncol = 1, nrow = eqm),
						matrix(0, ncol = n, nrow = eqm),
						matrix(0, ncol = n, nrow = eqm),
						matrix(0, ncol = n, nrow = eqm),
						matrix(0, ncol = n, nrow = eqm),
						matrix(0, ncol = n, nrow = eqm),
						matrix(0, ncol = n, nrow = eqm),
						matrix(0, ncol = n, nrow = eqm),
						matrix(0, ncol = 1, nrow = eqm),
						matrix(-optvars$eqB, ncol = 1, nrow = eqm) ) )
	} else{
		eqm = 0
	}
	
	# Constraint11: any user inequalities
	if(!is.null(optvars$ineq.mat)){
		ineqm = dim(optvars$ineq.mat)[1]
		Amat = rbind( Amat, cbind(
						optvars$ineq.mat,
						matrix(0, ncol = 1, nrow = ineqm),
						matrix(0, ncol = 1, nrow = ineqm),
						matrix(0, ncol = n, nrow = ineqm),
						matrix(0, ncol = n, nrow = ineqm),
						matrix(0, ncol = n, nrow = ineqm),
						matrix(0, ncol = n, nrow = ineqm),
						matrix(0, ncol = n, nrow = ineqm),
						matrix(0, ncol = n, nrow = ineqm),
						matrix(0, ncol = n, nrow = ineqm),
						matrix(0, ncol = 1, nrow = ineqm),
						matrix(-optvars$ineq.UB, ncol = 1, nrow = ineqm) ) )
		
		Amat = rbind( Amat, cbind(
						optvars$ineq.mat,
						matrix(0, ncol = 1, nrow = ineqm),
						matrix(0, ncol = 1, nrow = ineqm),
						matrix(0, ncol = n, nrow = ineqm),
						matrix(0, ncol = n, nrow = ineqm),
						matrix(0, ncol = n, nrow = ineqm),
						matrix(0, ncol = n, nrow = ineqm),
						matrix(0, ncol = n, nrow = ineqm),
						matrix(0, ncol = n, nrow = ineqm),
						matrix(0, ncol = n, nrow = ineqm),
						matrix(0, ncol = 1, nrow = ineqm),
						matrix(-optvars$ineq.LB, ncol = 1, nrow = ineqm) ) )
	} else{
		ineqm = 0
	}
	
	Amat = as.simple_triplet_matrix(Amat)
	
	rhs = c( 1, rep(0, n), rep(0, n), rep(0, 2*n), upa , 0, rep(0, m), rep(0, m), rep(0, eqm), rep(0, 2*ineqm) )	
	
	dir = c( "==", rep("==", 4*n), "==", rep("<=", 2 * ineqm), rep("==", eqm), "==", 
			rep("<=", m), rep(">=", m), rep("==", eqm), rep("<=", ineqm), rep(">=", ineqm) )
	
	# S, alpha, theta, u_t, v_t, phi_t, psi_t, z_t, slack1_t, slack2_t,  constant
	objL = c( rep(0, m), 1, 0, rep( probability/lpa, n ), rep(0, 6*n), 0 , 0 )
	
	xLB = c( rep(-Inf, m), 0, 0, rep(0, n), rep(0, n), rep(0, n), rep(0, n), rep(0, n), rep(0, 2*n), 1, 0.05)
	xUB = c( rep(Inf, m), Inf, Inf, rep(Inf, n), rep(Inf, n), rep(1, n), rep(1, 2 * n), rep(Inf, 2*n), 1, Inf)
	types = c(rep("C", m), "C", "C", rep("C", 4 * n), rep("B", n), rep("C", 2 * n), "C", "C")
	bounds = list( lower = list( ind = 1L:(4 + m + 7 * n), val = xLB ),
			upper = list( ind = 1L:(4 + m + 7 * n), val = xUB ) )
	
	ans = list(objL = objL, Amat = Amat, dir = dir, rhs = rhs, types = types, 
			bounds = bounds, dataidx = c(n, m), widx = 1:m)
	rm(Amat)
	rm(Data)
	gc()
	return(ans)
}


milplpm1.minrisk = function(optvars)
{
	# -[ Nxm ] -[ diag N ] + [1] = m+N+1
	Data = optvars$Data
	n = NROW(Data)
	m = NCOL(Data)
	colnames(Data) = NULL
	rownames(Data) = NULL
	probability = optvars$probability
	if(is.null(probability)) probability = rep(1/n, n)
	# -[ Nxm ] -[ diag N ] + [1] = m+N+1
	cons = lpmin.setup(optvars, cdimn = n+1)
	ineqcon = cons$ineqcon
	ineqB = cons$ineqB
	ineqn = cons$ineqn
	ineqm = cons$ineqm
	eqcon = cons$eqcon
	eqB = cons$eqB
	eqn = cons$eqn
	eqm = cons$eqm
	LB = optvars$LB
	UB = optvars$UB
	Amat = cbind( -Data, -diag(n), matrix(optvars$threshold, nrow = n, ncol = 1))
	Amat = rbind( Amat, ineqcon, eqcon )
	Amat = cbind( Amat, matrix(0, ncol = m, nrow = dim(Amat)[1]))
	Amat = rbind( Amat, cbind(-diag(m), matrix(0, ncol = n, nrow = m), matrix(0, ncol = 1, nrow = m), 
					diag(LB)) )
	Amat = rbind( Amat, cbind( diag(m), matrix(0, ncol = n, nrow = m), matrix(0, ncol = 1, nrow = m), 
					diag(-UB)) )
	Amat = rbind( Amat, cbind(matrix(0, ncol = m + n + 1, nrow = 1), matrix(1, ncol = m, nrow = 1) ) )	
	Amat = as.simple_triplet_matrix(Amat)
	
	rhs = rep(0, n)
	rhs = c( rhs, as.numeric(ineqB), as.numeric(eqB) , rep(0, 2*m), optvars$max.pos)
	dir = c( rep("<=", n), rep("<=", ineqn), rep("==", eqn), rep("<=", 2*m), "==")
	objL = c( rep(0, m), probability, 1 , rep(0, m))
	types = c( rep("C", m), rep("C", n), "C", rep("B", m))
	bounds = list( lower = list( ind = 1L:(m + n + 1 + m), val = c(LB,  rep(0, n), 1 , rep(0, m)) ),
			upper = list( ind = 1L:(m + n + 1 + m), val = c( UB, rep(Inf, n), 1, rep(1, m) ) ) )	
	ans = list(objL = objL, Amat = Amat, dir = dir, rhs = rhs, bounds = bounds, 
			problem.max = FALSE, types = types, 
			reward = optvars$mu, dataidx = c(n, m), widx = 1:m)
	rm(Amat)
	rm(Data)
	gc()
	return( ans )
}

milpmad.minrisk = function(optvars)
{
	Data = optvars$Data
	n = NROW(Data)
	m = NCOL(Data)
	colnames(Data) = NULL
	rownames(Data) = NULL
	probability = optvars$probability
	if(is.null(probability)) probability = rep(1/n, n)
	# -[ Nxm ] -[ diag N ] + [1] = m+N+1
	cons = lpmin.setup(optvars, cdimn = 2*n+1)
	ineqcon = cons$ineqcon
	ineqB = cons$ineqB
	ineqn = cons$ineqn
	ineqm = cons$ineqm
	eqcon = cons$eqcon
	eqB = cons$eqB
	eqn = cons$eqn
	eqm = cons$eqm
	LB = optvars$LB
	UB = optvars$UB
	
	Amat = cbind( -Data, matrix(optvars$benchmark, nrow = n, ncol = 1), diag(n), -diag(n) )
	Amat = rbind( Amat, ineqcon, eqcon )
	Amat = cbind( Amat, matrix(0, ncol = m, nrow = dim(Amat)[1]))

	Amat = rbind( Amat, cbind(-diag(m), matrix(0, ncol = 1, nrow = m), matrix(0, ncol = 2*n, nrow = m),
					diag(LB)) )
	Amat = rbind( Amat, cbind( diag(m), matrix(0, ncol = 1, nrow = m), matrix(0, ncol = 2*n, nrow = m),
					diag(-UB)) )
	Amat = rbind( Amat, cbind(matrix(0, ncol = m + 2*n + 1, nrow = 1), matrix(1, ncol = m, nrow = 1) ) )	
	Amat = as.simple_triplet_matrix(Amat)
	
	rhs = rep(0, n)
	
	rhs = c( rhs, as.numeric(ineqB), as.numeric(eqB) , rep(0, 2*m), optvars$max.pos)
	dir = c( rep("==", n), rep("<=", ineqn), rep("==", eqn) , rep("<=", 2*m), "==")
	objL = c( rep(0, m), 0, probability, probability, rep(0, m))
	bounds = list( lower = list( ind = 1L:(m + 2 * n + 1 + m), val = c(LB,  1, rep(0, 2 * n), rep(0, m) ) ),
			upper = list( ind = 1L:(m + 2 * n + 1 + m), val = c( UB, 1, rep(Inf, 2 * n), rep(1, m)) ) )	
	types = c( rep("C", m), "C", rep("C", 2*n), rep("B", m))
	
	ans = list(objL = objL, Amat = Amat, dir = dir, rhs = rhs, bounds = bounds, 
			problem.max = FALSE, reward = optvars$mu, 
			types = types, dataidx = c(n, m), widx = 1:m)
	rm(Amat)
	rm(Data)
	gc()
	return( ans )
}




milpminmax.minrisk = function(optvars)
{
	Data = optvars$Data
	n = NROW(Data)
	m = NCOL(Data)
	colnames(Data) = NULL
	rownames(Data) = NULL
	probability = optvars$probability
	if(is.null(probability)) probability = rep(1/n, n)
	# -[ Nxm ] -[ diag N ] + [1] = m+N+1
	cons = lpmin.setup(optvars, cdimn = 2)
	ineqcon = cons$ineqcon
	ineqB = cons$ineqB
	ineqn = cons$ineqn
	ineqm = cons$ineqm
	eqcon = cons$eqcon
	eqB = cons$eqB
	eqn = cons$eqn
	eqm = cons$eqm
	LB = optvars$LB
	UB = optvars$UB
	
	Amat = cbind( -Data, matrix(optvars$benchmark, nrow = n, ncol = 1), matrix(1, nrow = n, ncol = 1) )
	Amat = rbind( Amat, ineqcon, eqcon )
	Amat = cbind(Amat, matrix(0, ncol = m, nrow = dim(Amat)[1]))
	Amat = rbind( Amat, cbind(-diag(m), matrix(0, ncol = 2, nrow = m), diag(LB)) )
	Amat = rbind( Amat, cbind( diag(m), matrix(0, ncol = 2, nrow = m), diag(-UB)) )
	Amat = rbind( Amat, cbind(matrix(0, ncol = m + 2, nrow = 1), matrix(1, ncol = m, nrow = 1) ) )	
	Amat = as.simple_triplet_matrix(Amat)
	
	rhs = rep(0, n)
	rhs = c( rhs, as.numeric(ineqB), as.numeric(eqB), rep(0, 2*m), optvars$max.pos )
	
	dir = c( rep("<=", n), rep("<=", ineqn), rep("==", eqn), rep("<=", 2*m), "==")
	
	objL = c( rep(0, m), 0, -1 , rep(0, m))
	
	bounds = list( lower = list( ind = 1L:(m + 2 + m), val = c(LB,  1, -Inf, rep(0, m)) ),
			upper = list( ind = 1L:(m + 2 + m), val = c( UB, 1, Inf, rep(1, m)) ) )
	
	types = c( rep("C", m), rep("C", 2), rep("B", m))
	
	ans = list(objL = objL, Amat = Amat, dir = dir, rhs = rhs, bounds = bounds, 
			problem.max = FALSE, reward = optvars$mu,  types = types, 
			dataidx = c(n, m),  widx = 1:m, vidx = m+2)	
	rm(Amat)
	rm(Data)
	gc()
	return( ans )
}

milpcvar.minrisk = function(optvars)
{
	Data = optvars$Data
	n = NROW(Data)
	m = NCOL(Data)
	colnames(Data) = NULL
	rownames(Data) = NULL
	probability = optvars$probability
	if(is.null(probability)) probability = rep(1/n, n)
	# -[ Nxm ] -[ diag N ] + [1] = m+N+1
	cons = lpmin.setup(optvars, cdimn = n + 2)
	ineqcon = cons$ineqcon
	ineqB = cons$ineqB
	ineqn = cons$ineqn
	ineqm = cons$ineqm
	eqcon = cons$eqcon
	eqB = cons$eqB
	eqn = cons$eqn
	eqm = cons$eqm
	LB = optvars$LB
	UB = optvars$UB
	
	Amat = cbind( -Data,
			matrix(-1, nrow = n, ncol = 1), 
			-diag(n),
			matrix(optvars$benchmark, nrow = n, ncol = 1) )
	
	Amat = rbind( Amat, ineqcon, eqcon )
	Amat = cbind(Amat, matrix(0, ncol = m, nrow = dim(Amat)[1]))
	
	Amat = rbind( Amat, cbind(-diag(m), matrix(0, ncol = n+2, nrow = m), diag(LB)) )
	Amat = rbind( Amat, cbind( diag(m), matrix(0, ncol = n+2, nrow = m), diag(-UB)) )
	Amat = rbind( Amat, cbind(matrix(0, ncol = m + n + 2, nrow = 1), matrix(1, ncol = m, nrow = 1) ) )	
	Amat = as.simple_triplet_matrix(Amat)
	
	
	rhs = rep(0, n)
	rhs = c( rhs, as.numeric(ineqB), as.numeric(eqB), rep(0, 2*m), optvars$max.pos )
	
	dir = c( rep("<=", n), rep("<=", ineqn), rep("==", eqn), rep("<=", 2*m), "==")
	
	objL = c( rep(0, m), 1, probability/optvars$alpha, 0, rep(0, m))
	
	bounds = list( lower = list( ind = 1L:(m + n + 2 + m), val = c(LB,  -1, rep(0, n), 1, rep(0, m)) ),
			upper = list( ind = 1L:(m + n + 2 + m), val = c( UB, 1, rep(Inf, n), 1 , rep(1, m)) ) )
	
	types = c( rep("C", m), "C", rep("C", n), "C", rep("B", m))
	
	ans = list(objL = objL, Amat = Amat, dir = dir, rhs = rhs, bounds = bounds, 
			problem.max = FALSE, reward = optvars$mu, types = types, 
			dataidx = c(n, m), widx = 1:m, vidx = m+1)
	rm(Data)
	rm(Amat)
	gc()
	return( ans )	
}




milpcdar.minrisk = function(optvars)
{
	Data = optvars$Data
	n = NROW(Data)
	m = NCOL(Data)
	colnames(Data) = NULL
	rownames(Data) = NULL
	probability = optvars$probability
	if(is.null(probability)) probability = rep(1/n, n)
	# -[ Nxm ] -[ diag N ] + [1] = m+N+1
	cons = lpmin.setup(optvars, cdimn = 2*n + 2)
	ineqcon = cons$ineqcon
	ineqB = cons$ineqB
	ineqn = cons$ineqn
	ineqm = cons$ineqm
	eqcon = cons$eqcon
	eqB = cons$eqB
	eqn = cons$eqn
	eqm = cons$eqm
	LB = optvars$LB
	UB = optvars$UB
	
	xm = as.matrix( -diag(n) )
	idx = which(xm == -1, arr.ind = TRUE)
	xrow = idx[,1]
	xcol = idx[,2]
	xcol[2:length(xcol)] = xcol[2:length(xcol)] - 1
	diag(xm[ xrow[2:length(xrow)], xcol[2:length(xcol)] ]) = 1
	
	Amat = cbind( 
			-Data, 
			matrix(0, nrow = n, ncol = 1),
			matrix(0, nrow = n, ncol = n),
			as.matrix(xm),
			matrix(optvars$benchmark, nrow = n, ncol = 1))
	
	Amat = rbind( Amat, cbind(
					matrix(0, nrow = n, ncol = m), 
					matrix(-1, nrow = n, ncol = 1), 
					-1 * diag(n),
					diag(n),
					matrix(0, nrow = n, ncol = 1)) )
	
	Amat = rbind( Amat, ineqcon, eqcon )
	Amat = cbind(Amat, matrix(0, ncol = m, nrow = dim(Amat)[1]))
	
	Amat = rbind( Amat, cbind(-diag(m), matrix(0, ncol = 2*n+2, nrow = m), diag(LB)) )
	Amat = rbind( Amat, cbind( diag(m), matrix(0, ncol = 2*n+2, nrow = m), diag(-UB)) )
	Amat = rbind( Amat, cbind(matrix(0, ncol = m + 2*n + 2, nrow = 1), matrix(1, ncol = m, nrow = 1) ) )	
	Amat = as.simple_triplet_matrix(Amat)
	
	
	rhs = rep(0, 2 * n)
	rhs = c( rhs, as.numeric(ineqB), as.numeric(eqB), rep(0, 2*m), optvars$max.pos)
	
	dir = c( rep("<=", 2 * n), rep("<=", ineqn), rep("==", eqn), rep("<=", 2*m), "==")
	
	objL = c( rep(0, m), 1, probability/optvars$alpha, rep(0, n),  0, rep(0, m))
	
	bounds = list( lower = list( ind = 1L:(m + 2 * n + 2 + m), val = c( LB, -1, rep(0, 2 * n), 1, rep(0, m)) ),
			upper = list( ind = 1L:(m + 2 * n + 2 + m), val = c( UB, 1, rep(Inf, 2 * n), 1, rep(1, m) ) ) )	
	
	types = c( rep("C", m), "C", rep("C", 2*n), "C", rep("B", m))
	
	ans = list(objL = objL, Amat = Amat, dir = dir, rhs = rhs, bounds = bounds, 
			problem.max = FALSE, reward = optvars$mu, types = types, 
			dataidx = c(n, m), widx = 1:m, vidx = m+1)
	rm(Amat)
	rm(Data)
	gc()
	return( ans )	
}


milpport = function(optvars, ...)
{
	risk = c("mad", "minimax", "cvar", "cdar", "ev", "lpm", "lpmupm")[optvars$index[4]]
	setup = switch(risk,
				"mad"     = milpmad.minrisk(optvars),
				"minimax" = milpminmax.minrisk(optvars),
				"cvar"    = milpcvar.minrisk(optvars),
				"cdar"    = milpcdar.minrisk(optvars),
				"lpm"     = milplpm1.minrisk(optvars))
	sol = milpminsolver(setup, optvars, ...)
	
	#sol$risk = fun.risk(sol$weights, scenario, options, risk, benchmark)
	#sol$reward = sum(forecast * sol$weights)
	rm(setup)
	gc()
	return( sol )
}

milpminsolver = function(setup, optvars, ...)
{
	if( is.null(list(...)$verbose) ) verbose = FALSE else verbose = as.logical( list(...)$verbose )
	m = setup$dataidx[2]
	sol = try(Rglpk_solve_LP(obj = setup$objL, mat = setup$Amat, dir = setup$dir, 
					types = setup$types, rhs = setup$rhs, bounds = setup$bounds, 
					max = FALSE, verbose = verbose), silent = TRUE)
	if( inherits(sol, "try-error") ){
		status = "non-convergence"
		weights = rep(NA, m)
		risk = NA
		reward = NA
		if(optvars$index[4]==3) VaR = NA
		if(optvars$index[4]==4) DaR = NA
	} else{
		status  = sol$status
		weights = sol$solution[setup$widx]
		risk    = sol$optimum
		if( !is.null(setup$reward) ){
			reward = sum(setup$reward * weights)
		} else {
			reward = NULL
		}
		if(optvars$index[4]==3) VaR = sol$solution[setup$vidx]
		if(optvars$index[4]==4) DaR = sol$solution[setup$vidx]
	}
	rm(setup)
	gc()
	return( list(status = status, weights = weights, risk = risk, reward = reward, 
					multiplier = 1) )
}

# This is for the Rachev Ratio model
milpfracsolver = function(setup, optvars, ...)
{
	if( is.null(list(...)$verbose) ) verbose = FALSE else verbose = as.logical( list(...)$verbose )
	
	m = setup$dataidx[2]
	
	sol = try(Rglpk_solve_LP(obj = setup$objL, mat = setup$Amat, dir = setup$dir, 
					types = setup$type, rhs = setup$rhs, bounds = setup$bounds, 
					max = FALSE, verbose = verbose), silent = TRUE)
	ans = list()
	if( inherits(sol, "try-error") ){
		status = "non-convergence"
		weights = rep(NA, m)
		risk = NA
		reward = NA
		optimum = NA
		multiplier = NA
		VaR = NA
		DaR = NA
	} else{
		status = sol$status
		multiplier = sol$solution[setup$midx]
		tmp = sol$solution/multiplier
		weights = tmp[setup$widx]
		risk = sum(setup$objL * sol$solution)/multiplier
		reward = sum(optvars$mu * weights)
		optimum  = sol$optimum
		if(optvars$index[4]==3) VaR = sol$solution[setup$vidx]/multiplier else VaR = NA
		if(optvars$index[4]==4) DaR = sol$solution[setup$vidx]/multiplier else DaR = NA
	}
	
	
	
	rm(setup)
	
	return( list(status = status, weights = weights,  optimum = optimum) )
}

Try the parma package in your browser

Any scripts or data that you put into this service are public.

parma documentation built on Oct. 29, 2022, 1:08 a.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close