11.1 Portfolio Optimization¶
In this section the Markowitz portfolio optimization problem and variants are implemented using Optimizer API for Python.
Familiarity with Sec. 6.2 (From Linear to Conic Optimization) is recommended to follow the syntax used to create affine conic constraints (ACCs) throughout all the models appearing in this case study.
11.1.1 The Basic Model¶
The classical Markowitz portfolio optimization problem considers investing in
and covariance
The return of the investment is also a random variable
and variance
The standard deviation
is usually associated with risk.
The problem facing the investor is to rebalance the portfolio to achieve a good compromise between risk and expected return, e.g., maximize the expected return subject to a budget constraint and an upper bound (denoted
The variables
A popular choice is
Since
is the total investment. Clearly, the total amount invested must be equal to the initial wealth, which is
This leads to the first constraint
The second constraint
ensures that the variance, is bounded by the parameter
excludes the possibility of short-selling. This constraint can of course be excluded if short-selling is allowed.
The covariance matrix
In general the choice of
Hence, we may write the risk constraint as
or equivalently
where
Therefore, problem (11.1) can be written as
which is a conic quadratic optimization problem that can easily be formulated and solved with Optimizer API for Python. Subsequently we will use the example data
and
Using Cholesky factorization, this implies
In Sec. 11.1.3 (Factor model and efficiency), we present a different way of obtaining
Why a Conic Formulation?
Problem (11.1) is a convex quadratically constrained optimization problem that can be solved directly using MOSEK. Why then reformulate it as a conic quadratic optimization problem (11.3)? The main reason for choosing a conic model is that it is more robust and usually solves faster and more reliably. For instance it is not always easy to numerically validate that the matrix
Moreover, observe the constraint
more numerically robust than
for very small and very large values of
Example code
Listing 11.1 demonstrates how the basic Markowitz model (11.3) is implemented.
with mosek.Task() as task:
task.set_Stream(mosek.streamtype.log, sys.stdout.write)
# Holding variable x of length n
# No other auxiliary variables are needed in this formulation
task.appendvars(numvar)
# Optionally we can give the variables names
for j in range(0, n):
task.putvarname(voff_x + j, "x[%d]" % (1 + j))
# No short-selling in this model, all of x >= 0
task.putvarboundsliceconst(voff_x, n, mosek.boundkey.lo, 0.0, inf)
# One linear constraint: total budget
task.appendcons(1)
task.putconname(coff_bud, "budget")
task.putaijlist([coff_bud] * n, range(voff_x, voff_x + n), [1.0] * n) # e^T x
rtemp = w + sum(x0)
task.putconbound(coff_bud, mosek.boundkey.fx, rtemp, rtemp) # equals w + sum(x0)
# Input (gamma, GTx) in the AFE (affine expression) storage
# We need k+1 rows
task.appendafes(k + 1)
# The first affine expression = gamma
task.putafeg(0, gamma)
# The remaining k expressions comprise GT*x, we add them row by row
# In more realisic scenarios it would be better to extract nonzeros and input in sparse form
for i in range(0, k):
task.putafefrow(i + 1, range(voff_x, voff_x + n), GT[i])
# Input the affine conic constraint (gamma, GT*x) \in QCone
# Add the quadratic domain of dimension k+1
qdom = task.appendquadraticconedomain(k + 1)
# Add the constraint
task.appendaccseq(qdom, 0, None)
task.putaccname(0, "risk")
# Objective: maximize expected return mu^T x
task.putclist(range(voff_x, voff_x + n), mu)
task.putobjsense(mosek.objsense.maximize)
# Dump the problem to a human readable PTF file.
task.writedata("dump.ptf")
# Solve the problem
task.optimize()
# Display solution summary for quick inspection of results.
# In this simplified example we skip checks for problem and solution status
task.solutionsummary(mosek.streamtype.msg)
# Check if the interior point solution is an optimal point
solsta = task.getsolsta(mosek.soltype.itr)
if (solsta != mosek.solsta.optimal):
# See https://docs.mosek.com/latest/pythonapi/accessing-solution.html about handling solution statuses.
raise Exception(f"Unexpected solution status: {solsta}")
# Retrieve results
xx = task.getxxslice(mosek.soltype.itr, voff_x, voff_x + n)
expret = task.getprimalobj(mosek.soltype.itr)
print(f'Expected return: {expret:.10e} Std. deviation: {gamma:.4e}')
np.set_printoptions(precision=4)
print(f'Optimal portfolio: {np.array(xx)}')
The code is organized as follows:
We have
optimization variables, one per each asset in the portfolio. They correspond to the variable from (11.1) and their indices as variables in the task are from to (inclusive).The linear part of the problem: budget constraint, no-short-selling bounds and the objective are added in the linear data of the task (
matrix, vector and bounds) following the techniques introduced in the tutorial of Sec. 6.1 (Linear Optimization).For the quadratic constraint we follow the path introduced in the tutorial of Sec. 6.2 (From Linear to Conic Optimization). We add the vector
to the affine expression storage (AFE), create a quadratic domain of suitable length, and add the affine conic constraint (ACC) with the selected affine expressions. In the segment# Input the affine conic constraint (gamma, GT*x) \in QCone # Add the quadratic domain of dimension k+1 qdom = task.appendquadraticconedomain(k + 1) # Add the constraint task.appendaccseq(qdom, 0, None)
we use
Task.appendaccseq
to append a single ACC with the quadratic domainqdom
and with a sequence of affine expressions starting at position in the AFE storage and of length equal to the dimension ofqdom
. This is the simplest way to achieve what we need, since previously we also stored the required rows in AFE in the same order.
11.1.2 The Efficient Frontier¶
The portfolio computed by the Markowitz model is efficient in the sense that there is no other portfolio giving a strictly higher return for the same amount of risk. An efficient portfolio is also sometimes called a Pareto optimal portfolio. Clearly, an investor should only invest in efficient portfolios and therefore it may be relevant to present the investor with all efficient portfolios so the investor can choose the portfolio that has the desired tradeoff between return and risk.
Given a nonnegative
is one standard way to trade the expected return against penalizing variance. Note that, in contrast to the previous example, we explicitly use the variance (
The parameter

Fig. 11.1 The efficient frontier for the sample data.¶
Example code
Listing 11.2 demonstrates how to compute the efficient portfolios for several values of
here
to download.¶ with mosek.Task() as task:
task.set_Stream(mosek.streamtype.log, sys.stdout.write)
# Variables:
task.appendvars(numvar)
# Optionally we can give the variables names
for j in range(0, n):
task.putvarname(voff_x + j, "x[%d]" % (1 + j))
task.putvarname(voff_s, "s")
# No short-selling in this model, all of x >= 0
task.putvarboundsliceconst(voff_x, n, mosek.boundkey.lo, 0.0, inf)
# s is free variable
task.putvarbound(voff_s, mosek.boundkey.fr, -inf, inf)
# One linear constraint: total budget
task.appendcons(1)
task.putconname(coff_bud, "budget")
task.putaijlist([coff_bud] * n, range(voff_x, voff_x + n), [1.0] * n) # e^T x
rtemp = w + sum(x0)
task.putconbound(coff_bud, mosek.boundkey.fx, rtemp, rtemp) # equals w + sum(x0)
# Input (gamma, GTx) in the AFE (affine expression) storage
# We build the following F and g for variables [x, s]:
# [0, 1] [0 ]
# F = [0, 0], g = [0.5]
# [GT,0] [0 ]
# We need k+2 rows
task.appendafes(k + 2)
# The first affine expression is variable s (last variable, index n)
task.putafefrow(0, [n], [1.0])
# The second affine expression is constant 0.5
task.putafeg(1, 0.5)
# The remaining k affine expressions comprise GT*x, we add them row by row
# In more realisic scenarios it would be better to extract nonzeros and input in sparse form
for i in range(0, k):
task.putafefrow(i + 2, range(voff_x, voff_x + n), GT[i])
# Input the affine conic constraint (s, 0.5, GT*x) \in RQCone
# Add the quadratic domain of dimension k+1
rqdom = task.appendrquadraticconedomain(k + 2)
# Add the constraint
task.appendaccseq(rqdom, 0, None)
task.putaccname(0, "risk")
# Set objective coefficients (x part): mu'x - alpha * s
task.putclist(range(voff_x, voff_x + n), mu)
task.putobjsense(mosek.objsense.maximize)
# Turn all log output off.
task.putintparam(mosek.iparam.log, 0)
for alpha in alphas:
# Dump the problem to a human readable PTF file.
task.writedata("dump.ptf")
task.putcj(voff_s, -alpha)
task.optimize()
# Display the solution summary for quick inspection of results.
# task.solutionsummary(mosek.streamtype.msg)
# Check if the interior point solution is an optimal point
solsta = task.getsolsta(mosek.soltype.itr)
if solsta in [mosek.solsta.optimal]:
expret = 0.0
x = task.getxxslice(mosek.soltype.itr, voff_x, voff_x + n)
for j in range(0, n):
expret += mu[j] * x[j]
stddev = np.sqrt(task.getxxslice(mosek.soltype.itr, voff_s, voff_s + 1))
print("alpha = {0:.2e} exp. ret. = {1:.3e}, std. dev. {2:.3e}".format(alpha, expret, stddev[0]))
else:
# See https://docs.mosek.com/latest/pythonapi/accessing-solution.html about handling solution statuses.
print("An error occurred when solving for alpha=%e\n" % alpha)
raise Exception(f"Unexpected solution status: {solsta}")
Note that we changed the coefficient
11.1.3 Factor model and efficiency¶
In practice it is often important to solve the portfolio problem very quickly. Therefore, in this section we discuss how to improve computational efficiency at the modeling stage.
The computational cost is of course to some extent dependent on the number of constraints and variables in the optimization problem. However, in practice a more important factor is the sparsity: the number of nonzeros used to represent the problem. Indeed it is often better to focus on the number of nonzeros in
In other words if the computational efficiency should be improved then it is always good idea to start with focusing at the covariance matrix. As an example assume that
where
One possible choice for
because then
This choice requires storage proportional to
The example above exploits the so-called factor structure and demonstrates that an alternative choice of
The lesson to be learned is that it is important to investigate how the covariance matrix is formed. Given this knowledge it might be possible to make a special choice for
Factor model in finance
Factor model structure is typical in financial context. It is common to model security returns as the sum of two components using a factor model. The first component is the linear combination of a small number of factors common among a group of securities. The second component is a residual, specific to each security. It can be written as
Such a model will result in the covariance structure
where
Example code
Here we will work with the example data of a two-factor model (
and the factor covariance matrix is
giving
Then the matrix
This matrix is indeed very sparse.
In general, we get an
Example code
In the following we demonstrate how to write code to compute the matrix
B = np.array([
[0.4256, 0.1869],
[0.2413, 0.3877],
[0.2235, 0.3697],
[0.1503, 0.4612],
[1.5325, -0.2633],
[1.2741, -0.2613],
[0.6939, 0.2372],
[0.5425, 0.2116]
])
S_F = np.array([
[0.0620, 0.0577],
[0.0577, 0.0908]
])
theta = np.array([0.0720, 0.0508, 0.0377, 0.0394, 0.0663, 0.0224, 0.0417, 0.0459])
Then the matrix
P = np.linalg.cholesky(S_F)
G_factor = B @ P
The code for computing an optimal portfolio in the factor model is very similar to the one from the basic model in Listing 11.1 with one notable exception: we construct the expression
# Input (gamma, G_factor_T x, diag(sqrt(theta))*x) in the AFE (affine expression) storage
# We need k+n+1 rows and we fill them in in three parts
task.appendafes(k + n + 1)
# 1. The first affine expression = gamma, will be specified later
# 2. The next k expressions comprise G_factor_T*x, we add them row by row
# transposing the matrix G_factor on the fly
for i in range(0, k):
task.putafefrow(i + 1, range(voff_x, voff_x + n), np.array(G_factor[:,i]))
# 3. The remaining n rows contain sqrt(theta) on the diagonal
task.putafefentrylist(range(k + 1, k + 1 + n), range(voff_x, voff_x + n), np.sqrt(theta))
The full code is demonstrated below:
import mosek
import sys
import numpy as np
if __name__ == '__main__':
n = 8
w = 1.0
mu = [0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379]
x0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
B = np.array([
[0.4256, 0.1869],
[0.2413, 0.3877],
[0.2235, 0.3697],
[0.1503, 0.4612],
[1.5325, -0.2633],
[1.2741, -0.2613],
[0.6939, 0.2372],
[0.5425, 0.2116]
])
S_F = np.array([
[0.0620, 0.0577],
[0.0577, 0.0908]
])
theta = np.array([0.0720, 0.0508, 0.0377, 0.0394, 0.0663, 0.0224, 0.0417, 0.0459])
P = np.linalg.cholesky(S_F)
G_factor = B @ P
k = G_factor.shape[1]
gammas = [0.24, 0.28, 0.32, 0.36, 0.4, 0.44, 0.48]
inf = 0.0 # This value has no significance
# Variable offsets
numvar = n
voff_x = 0
# Constraints offsets
numcon = 1
coff_bud = 0
with mosek.Task() as task:
task.set_Stream(mosek.streamtype.log, sys.stdout.write)
# Holding variable x of length n
# No other auxiliary variables are needed in this formulation
task.appendvars(numvar)
# Optionally we can give the variables names
for j in range(0, n):
task.putvarname(voff_x + j, "x[%d]" % (1 + j))
# No short-selling in this model, all of x >= 0
task.putvarboundsliceconst(voff_x, n, mosek.boundkey.lo, 0.0, inf)
# One linear constraint: total budget
task.appendcons(1)
task.putconname(coff_bud, "budget")
task.putaijlist([coff_bud] * n, range(voff_x, voff_x + n), [1.0] * n) # e^T x
rtemp = w + sum(x0)
task.putconbound(coff_bud, mosek.boundkey.fx, rtemp, rtemp) # equals w + sum(x0)
# Input (gamma, G_factor_T x, diag(sqrt(theta))*x) in the AFE (affine expression) storage
# We need k+n+1 rows and we fill them in in three parts
task.appendafes(k + n + 1)
# 1. The first affine expression = gamma, will be specified later
# 2. The next k expressions comprise G_factor_T*x, we add them row by row
# transposing the matrix G_factor on the fly
for i in range(0, k):
task.putafefrow(i + 1, range(voff_x, voff_x + n), np.array(G_factor[:,i]))
# 3. The remaining n rows contain sqrt(theta) on the diagonal
task.putafefentrylist(range(k + 1, k + 1 + n), range(voff_x, voff_x + n), np.sqrt(theta))
# Input the affine conic constraint (gamma, G_factor_T x, diag(sqrt(theta))*x) \in QCone
# Add the quadratic domain of dimension k+n+1
qdom = task.appendquadraticconedomain(k + n + 1)
# Add the constraint
task.appendaccseq(qdom, 0, None)
task.putaccname(0, "risk")
# Objective: maximize expected return mu^T x
task.putclist(range(voff_x, voff_x + n), mu)
task.putobjsense(mosek.objsense.maximize)
for gamma in gammas:
# Specify gamma in ACC
task.putafeg(0, gamma)
# Dump the problem to a human readable PTF file.
task.writedata("dump.ptf")
# Solve the problem
task.optimize()
# Display solution summary for quick inspection of results.
# In this simplified example we skip checks for problem and solution status
task.solutionsummary(mosek.streamtype.msg)
# Check if the interior point solution is an optimal point
solsta = task.getsolsta(mosek.soltype.itr)
if (solsta != mosek.solsta.optimal):
# See https://docs.mosek.com/latest/pythonapi/accessing-solution.html about handling solution statuses.
raise Exception(f"Unexpected solution status: {solsta}")
# Retrieve results
xx = task.getxxslice(mosek.soltype.itr, voff_x, voff_x + n)
expret = task.getprimalobj(mosek.soltype.itr)
print(f'Expected return: {expret:.10e} Std. deviation: {gamma:.4e}')
np.set_printoptions(precision=4)
print(f'Optimal portfolio: {np.array(xx)}')
11.1.4 Slippage Cost¶
The basic Markowitz model assumes that there are no costs associated with trading the assets and that the returns of the assets are independent of the amount traded. Neither of those assumptions is usually valid in practice. Therefore, a more realistic model is
Here
and
11.1.5 Market Impact Costs¶
If the initial wealth is fairly small and no short selling is allowed, then the holdings will be small and the traded amount of each asset must also be small. Therefore, it is reasonable to assume that the prices of the assets are independent of the amount traded. However, if a large volume of an asset is sold or purchased, the price, and hence return, can be expected to change. This effect is called market impact costs. It is common to assume that the market impact cost for asset
where
Hence, it follows that
Unfortunately this set of constraints is nonconvex due to the constraint
but in many cases the constraint may be replaced by the relaxed constraint
For instance if the universe of assets contains a risk free asset then
cannot hold for an optimal solution.
If the optimal solution has the property (11.9) then the market impact cost within the model is larger than the true market impact cost and hence money are essentially considered garbage and removed by generating transaction costs. This may happen if a portfolio with very small risk is requested because the only way to obtain a small risk is to get rid of some of the assets by generating transaction costs. We generally assume that this is not the case and hence the models (11.7) and (11.8) are equivalent.
The above observations lead to
The revised budget constraint
specifies that the initial wealth covers the investment and the transaction costs. It should be mentioned that transaction costs of the form
where
See the Modeling Cookbook for details.
Example code
Listing 11.5 demonstrates how to compute an optimal portfolio when market impact cost are included.
with mosek.Task() as task:
task.set_Stream(mosek.streamtype.log, sys.stdout.write)
# Variables (vector of x, c, z)
task.appendvars(numvar)
for j in range(0, n):
task.putvarname(voff_x + j, "x[%d]" % (j + 1))
task.putvarname(voff_c + j, "c[%d]" % (j + 1))
task.putvarname(voff_z + j, "z[%d]" % (j + 1))
# Apply variable bounds (x >= 0, c and z free)
task.putvarboundsliceconst(voff_x, voff_x + n, mosek.boundkey.lo, 0.0, inf)
task.putvarboundsliceconst(voff_c, voff_c + n, mosek.boundkey.fr, -inf, inf)
task.putvarboundsliceconst(voff_z, voff_z + n, mosek.boundkey.fr, -inf, inf)
# Linear constraints
# - Budget
task.appendcons(1)
task.putconname(coff_bud, "budget")
task.putaijlist([coff_bud] * n, range(voff_x, voff_x + n), [1.0] * n) # e^T x
task.putaijlist([coff_bud] * n, range(voff_c, voff_c + n), m) # m^T c
rtemp = w + sum(x0)
task.putconbound(coff_bud, mosek.boundkey.fx, rtemp, rtemp) # equals w + sum(x0)
# - Absolute value
task.appendcons(2 * n)
for i in range(0, n):
task.putconname(coff_abs1 + i, "zabs1[%d]" % (1 + i))
task.putconname(coff_abs2 + i, "zabs2[%d]" % (1 + i))
task.putaijlist(range(coff_abs1, coff_abs1 + n), range(voff_x, voff_x + n), n * [-1.0])
task.putaijlist(range(coff_abs1, coff_abs1 + n), range(voff_z, voff_z + n), n * [1.0])
task.putconboundlist(range(coff_abs1, coff_abs1 + n), [mosek.boundkey.lo] * n, [-x0[j] for j in range(0, n)], [inf] * n)
task.putaijlist(range(coff_abs2, coff_abs2 + n), range(voff_x, voff_x + n), n * [1.0])
task.putaijlist(range(coff_abs2, coff_abs2 + n), range(voff_z, voff_z + n), n * [1.0])
task.putconboundlist(range(coff_abs2, coff_abs2 + n), [mosek.boundkey.lo] * n, x0, [inf] * n)
# ACCs
aoff_q = 0
aoff_pow = k + 1
# - (gamma, GTx) in Q(k+1)
# The part of F and g for variable x:
# [0, 0, 0] [gamma]
# F = [GT, 0, 0], g = [0 ]
task.appendafes(k + 1)
task.putafeg(aoff_q, gamma)
for i in range(0, k):
task.putafefrow(aoff_q + 1 + i, range(voff_x, voff_x + n), GT[i])
qdom = task.appendquadraticconedomain(k + 1)
task.appendaccseq(qdom, aoff_q, None)
task.putaccname(aoff_q, "risk")
# - (c_j, 1, z_j) in P3(2/3, 1/3)
# The part of F and g for variables [c, z]:
# [0, I, 0] [0]
# F = [0, 0, I], g = [0]
# [0, 0, 0] [1]
task.appendafes(2 * n + 1)
task.putafefentrylist(range(aoff_pow, aoff_pow + n), range(voff_c, voff_c + n), [1.0] * n)
task.putafefentrylist(range(aoff_pow + n, aoff_pow + 2 * n), range(voff_z, voff_z + n), [1.0] * n)
task.putafeg(aoff_pow + 2 * n, 1.0)
# We use one row from F and g for both c_j and z_j, and the last row of F and g for the constant 1.
# NOTE: Here we reuse the last AFE and the power cone n times, but we store them only once.
powdom = task.appendprimalpowerconedomain(3, [2, 1])
afe_list = [(aoff_pow + i, aoff_pow + 2 * n, aoff_pow + n + i) for i in range(0, n)]
flat_afe_list = [idx for sublist in afe_list for idx in sublist]
task.appendaccs([powdom] * n, flat_afe_list, None)
for i in range(0, n):
task.putaccname(i + 1, "market_impact[%d]" % i)
# Objective
task.putclist(range(voff_x, voff_x + n), mu)
task.putobjsense(mosek.objsense.maximize)
# Turn all log output off.
# task.putintparam(mosek.iparam.log,0)
# Dump the problem to a human readable PTF file.
task.writedata("dump.ptf")
task.optimize()
# Display the solution summary for quick inspection of results.
task.solutionsummary(mosek.streamtype.msg)
# Check if the interior point solution is an optimal point
solsta = task.getsolsta(mosek.soltype.itr)
if (solsta != mosek.solsta.optimal):
# See https://docs.mosek.com/latest/pythonapi/accessing-solution.html about handling solution statuses.
raise Exception(f"Unexpected solution status: {solsta}")
expret = 0.0
x = task.getxxslice(mosek.soltype.itr, voff_x, voff_x + n)
for j in range(0, n):
expret += mu[j] * x[j]
print("\nExpected return %e for gamma %e\n" % (expret, gamma))
Note that in the following part of the code:
task.putafeg(aoff_pow + 2 * n, 1.0)
# We use one row from F and g for both c_j and z_j, and the last row of F and g for the constant 1.
# NOTE: Here we reuse the last AFE and the power cone n times, but we store them only once.
powdom = task.appendprimalpowerconedomain(3, [2, 1])
afe_list = [(aoff_pow + i, aoff_pow + 2 * n, aoff_pow + n + i) for i in range(0, n)]
flat_afe_list = [idx for sublist in afe_list for idx in sublist]
task.appendaccs([powdom] * n, flat_afe_list, None)
for i in range(0, n):
task.putaccname(i + 1, "market_impact[%d]" % i)
we create a sequence of power cones of the form
Moreover, note that the second coordinate of all these affine conic constraints is the same affine expression equal to aoff_pow + 2 * n
) and reuse it in all the ACCs.
11.1.6 Transaction Costs¶
Now assume there is a cost associated with trading asset
Hence, whenever asset
(11.11)¶
First observe that
We choose
Example code
The following example code demonstrates how to compute an optimal portfolio when transaction costs are included.
with mosek.Task() as task:
task.set_Stream(mosek.streamtype.log, sys.stdout.write)
# Variables (vector of x, z, y)
task.appendvars(numvar)
for j in range(0, n):
task.putvarname(voff_x + j, "x[%d]" % (j + 1))
task.putvarname(voff_z + j, "z[%d]" % (j + 1))
task.putvarname(voff_y + j, "y[%d]" % (j + 1))
# Apply variable bounds (x >= 0, z free, y binary)
task.putvarboundsliceconst(voff_x, voff_x + n, mosek.boundkey.lo, 0.0, inf)
task.putvarboundsliceconst(voff_z, voff_z + n, mosek.boundkey.fr, -inf, inf)
task.putvarboundsliceconst(voff_y, voff_y + n, mosek.boundkey.ra, 0.0, 1.0)
task.putvartypelist(range(voff_y, voff_y + n), [mosek.variabletype.type_int] * n)
# Linear constraints
# - Budget
task.appendcons(1)
task.putconname(coff_bud, "budget")
task.putaijlist([0] * n, range(voff_x, voff_x + n), [1.0] * n) # e^T x
task.putaijlist([0] * n, range(voff_z, voff_z + n), g) # g^T z
task.putaijlist([0] * n, range(voff_y, voff_y + n), f) # f^T y
U = w + sum(x0)
task.putconbound(coff_bud, mosek.boundkey.fx, U, U) # equals w + sum(x0)
# - Absolute value
task.appendcons(2 * n)
for i in range(0, n):
task.putconname(coff_abs1 + i, "zabs1[%d]" % (1 + i))
task.putconname(coff_abs2 + i, "zabs2[%d]" % (1 + i))
task.putaijlist(range(coff_abs1, coff_abs1 + n), range(voff_x, voff_x + n), [-1.0] * n)
task.putaijlist(range(coff_abs1, coff_abs1 + n), range(voff_z, voff_z + n), [1.0] * n)
task.putconboundlist(range(coff_abs1, coff_abs1 + n), [mosek.boundkey.lo] * n, [-x0[j] for j in range(0, n)], [inf] * n)
task.putaijlist(range(coff_abs2, coff_abs2 + n), range(voff_x, voff_x + n), [1.0] * n)
task.putaijlist(range(coff_abs2, coff_abs2 + n), range(voff_z, voff_z + n), [1.0] * n)
task.putconboundlist(range(coff_abs2, coff_abs2 + n), [mosek.boundkey.lo] * n, x0, [inf] * n)
# - Switch
task.appendcons(n)
for i in range(0, n):
task.putconname(coff_swi + i, "switch[%d]" % (1 + i))
task.putaijlist(range(coff_swi, coff_swi + n), range(voff_z, voff_z + n), [1.0] * n)
task.putaijlist(range(coff_swi, coff_swi + n), range(voff_y, voff_y + n), [-U] * n)
task.putconboundlist(range(coff_swi, coff_swi + n), [mosek.boundkey.up] * n, [-inf] * n, [0.0] * n)
# ACCs
aoff_q = 0
# - (gamma, GTx) in Q(k+1)
# The part of F and g for variable x:
# [0, 0, 0] [gamma]
# F = [GT, 0, 0], g = [0 ]
task.appendafes(k + 1)
task.putafeg(aoff_q, gamma)
for i in range(0, k):
task.putafefrow(aoff_q + i + 1, range(voff_x, voff_x + n), GT[i])
qdom = task.appendquadraticconedomain(k + 1)
task.appendaccseq(qdom, aoff_q, None)
task.putaccname(0, "risk")
# Objective
task.putclist(range(voff_x, voff_x + n), mu)
task.putobjsense(mosek.objsense.maximize)
# Turn all log output off.
# task.putintparam(mosek.iparam.log,0)
# Dump the problem to a human readable PTF file.
task.writedata("dump.ptf")
task.optimize()
# Display the solution summary for quick inspection of results.
task.solutionsummary(mosek.streamtype.msg)
# Check if the integer solution is an optimal point
solsta = task.getsolsta(mosek.soltype.itg)
if (solsta != mosek.solsta.integer_optimal):
# See https://docs.mosek.com/latest/pythonapi/accessing-solution.html about handling solution statuses.
raise Exception(f"Unexpected solution status: {solsta}")
expret = 0.0
x = task.getxxslice(mosek.soltype.itg, voff_x, voff_x + n)
for j in range(0, n):
expret += mu[j] * x[j]
tcost = 0.0
z = task.getxxslice(mosek.soltype.itg, voff_z, voff_z + n)
y = task.getxxslice(mosek.soltype.itg, voff_y, voff_y + n)
for j in range(0, n):
tcost += g[j] * z[j] + f[j] * y[j]
print("\nExpected return %e for gamma %e. Transaction cost: %e\n" % (expret, gamma, tcost))
11.1.7 Cardinality constraints¶
Another method to reduce costs involved with processing transactions is to only change positions in a small number of assets. In other words, at most
This type of constraint can be again modeled by introducing a binary variable
(11.12)¶
were
Example code
The following example code demonstrates how to compute an optimal portfolio with cardinality bounds.
def markowitz_with_card(n, k, x0, w, gamma, mu, GT, K):
with mosek.Task() as task:
task.set_Stream(mosek.streamtype.log, sys.stdout.write)
# Offset of variables.
numvar = 3 * n
voff_x, voff_z, voff_y = 0, n, 2 * n
# Offset of constraints.
numcon = 3 * n + 2
coff_bud, coff_abs1, coff_abs2, coff_swi, coff_card = 0, 1, 1 + n, 1 + 2 * n, 1 + 3 * n
# Variables (vector of x, z, y)
task.appendvars(numvar)
for j in range(0, n):
task.putvarname(voff_x + j, "x[%d]" % (j + 1))
task.putvarname(voff_z + j, "z[%d]" % (j + 1))
task.putvarname(voff_y + j, "y[%d]" % (j + 1))
# Apply variable bounds (x >= 0, z free, y binary)
task.putvarboundsliceconst(voff_x, voff_x + n, mosek.boundkey.lo, 0.0, inf)
task.putvarboundsliceconst(voff_z, voff_z + n, mosek.boundkey.fr, -inf, inf)
task.putvarboundsliceconst(voff_y, voff_y + n, mosek.boundkey.ra, 0.0, 1.0)
task.putvartypelist(range(voff_y, voff_y + n), [mosek.variabletype.type_int] * n)
# Linear constraints
# - Budget
task.appendcons(1)
task.putconname(coff_bud, "budget")
task.putaijlist([coff_bud] * n, range(voff_x, voff_x + n), [1.0] * n) # e^T x
U = w + sum(x0)
task.putconbound(coff_bud, mosek.boundkey.fx, U, U) # = w + sum(x0)
# - Absolute value
task.appendcons(2 * n)
for i in range(0, n):
task.putconname(coff_abs1 + i, "zabs1[%d]" % (1 + i))
task.putconname(coff_abs2 + i, "zabs2[%d]" % (1 + i))
task.putaijlist(range(coff_abs1, coff_abs1 + n), range(voff_x, voff_x + n), [-1.0] * n)
task.putaijlist(range(coff_abs1, coff_abs1 + n), range(voff_z, voff_z + n), [1.0] * n)
task.putconboundlist(range(coff_abs1, coff_abs1 + n), [mosek.boundkey.lo] * n, [-x0[j] for j in range(0, n)], [inf] * n)
task.putaijlist(range(coff_abs2, coff_abs2 + n), range(voff_x, voff_x + n), [1.0] * n)
task.putaijlist(range(coff_abs2, coff_abs2 + n), range(voff_z, voff_z + n), [1.0] * n)
task.putconboundlist(range(coff_abs2, coff_abs2 + n), [mosek.boundkey.lo] * n, x0, [inf] * n)
# - Switch
task.appendcons(n)
for i in range(0, n):
task.putconname(coff_swi + i, "switch[%d]" % (1 + i))
task.putaijlist(range(coff_swi, coff_swi + n), range(voff_z, voff_z + n), [1.0] * n)
task.putaijlist(range(coff_swi, coff_swi + n), range(voff_y, voff_y + n), [-U] * n)
task.putconboundlist(range(coff_swi, coff_swi + n), [mosek.boundkey.up] * n, [-inf] * n, [0.0] * n)
# - Cardinality
task.appendcons(1)
task.putconname(coff_card, "cardinality")
task.putaijlist([coff_card] * n, range(voff_y, voff_y + n), [1.0] * n) # e^T y
task.putconbound(coff_card, mosek.boundkey.up, -inf, K) # <= K
# ACCs
aoff_q = 0
# - (gamma, GTx) in Q(k+1)
# The part of F and g for variable x:
# [0, 0, 0] [gamma]
# F = [GT, 0, 0], g = [0 ]
task.appendafes(k + 1)
task.putafeg(aoff_q, gamma)
for i in range(0, k):
task.putafefrow(aoff_q + i + 1, range(voff_x, voff_x + n), GT[i])
qdom = task.appendquadraticconedomain(k + 1)
task.appendaccseq(qdom, aoff_q, None)
task.putaccname(0, "risk")
# Objective
task.putclist(range(voff_x, voff_x + n), mu)
task.putobjsense(mosek.objsense.maximize)
# Turn all log output off.
task.putintparam(mosek.iparam.log,0)
# Dump the problem to a human readable PTF file.
task.writedata("dump.ptf")
task.optimize()
# Check if the integer solution is an optimal point
solsta = task.getsolsta(mosek.soltype.itg)
if (solsta != mosek.solsta.integer_optimal):
# See https://docs.mosek.com/latest/pythonapi/accessing-solution.html about handling solution statuses.
raise Exception(f"Unexpected solution status: {solsta}")
# Display the solution summary for quick inspection of results.
#task.solutionsummary(mosek.streamtype.msg)
return task.getxxslice(mosek.soltype.itg, voff_x + 0, voff_x + n)
If we solve our running example with
Bound 1 Solution: 0.0000e+00 0.0000e+00 1.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00
Bound 2 Solution: 0.0000e+00 0.0000e+00 3.5691e-01 0.0000e+00 0.0000e+00 6.4309e-01 -0.0000e+00 0.0000e+00
Bound 3 Solution: 0.0000e+00 0.0000e+00 1.9258e-01 0.0000e+00 0.0000e+00 5.4592e-01 2.6150e-01 0.0000e+00
Bound 4 Solution: 0.0000e+00 0.0000e+00 2.0391e-01 0.0000e+00 6.7098e-02 4.9181e-01 2.3718e-01 0.0000e+00
Bound 5 Solution: 0.0000e+00 3.1970e-02 1.7028e-01 0.0000e+00 7.0741e-02 4.9551e-01 2.3150e-01 0.0000e+00
Bound 6 Solution: 0.0000e+00 3.1970e-02 1.7028e-01 0.0000e+00 7.0740e-02 4.9551e-01 2.3150e-01 0.0000e+00
Bound 7 Solution: 0.0000e+00 3.1970e-02 1.7028e-01 0.0000e+00 7.0740e-02 4.9551e-01 2.3150e-01 0.0000e+00
Bound 8 Solution: 1.9557e-10 2.6992e-02 1.6706e-01 2.9676e-10 7.1245e-02 4.9559e-01 2.2943e-01 9.6905e-03