cvxopt optimization

Hi, wondering if anyone might have an idea how to handle this case when it comes to portfolio optimization and how to construct the matrices. Simple case is fine using portfolios = [solvers.qp(mu*P, -q, G, h, A, b)['x'] for mu in mus] where sum(x)=1 and x>=0 but what happens if you have n assets and you set upper and lower bounds. Example, if I have 3 assets and say lower bound is 40% and upper bound is 100%, the optimizer will fail because lower weights on all 3 assets will sum to greater than 1. Therefore, that means, 1 asset must fall out. Only way I can understand this case is optimize on the combinations of 3 assets with same lower bound of 40%. Thus, optimize(0,1) then (0,2) and finally (1,2) and find the combination with lowest stdev if that is your objective but if assets are large this process will be computationally heavy. Each combination can have lower bound of 40% and sum to 1. So, my question, how do you change G, h, A, b so this will work. I think you need to change lower bound <= x <= upper bound to an equality but can not figure this out. Any feedback greatly appreciated. thanks

13 responses

I don't understand? One is an equality constraint (sum all to one) and another is an inequality constraint (bounds). Where is the issue? use A and b for equality and G and h for inequality.

here is a quick example using cvxpy which is similar to cvxopt
n=3 (3 assets)
prob = Problem(Maximize(ret - gamma*risk), constraints)
constraints=[sum_entries(w) == 1, w >= 0.30, w<=1.0]

weights
[matrix([[ 0.4], [ 0.3],
[ 0.3]])]
no problem

but....
constraints=[sum_entries(w) == 1, w >= 0.40, w<=1.0]
weights=[None]

with three assets and say lower bound is 40%, solution can not be found. Not sure how to alter the matrices

When working with optimizers and their matrices, I find I always have to write out the constraints I want mathematically, and then write out the equation matrices, and only then try to code it. Anything else and I confuse myself. I use a http://www.amazon.com/uni-ball-KuruToga-Mechanical-Starter-1751934/dp/B0026ICM1E for this.

constraints=[sum_entries(w) == 1, w >= 0.40, w<=1.0]
to
constraints=[sum_entries(w) <= 1, w >= 0.40, w<=1.0]

and if the answer sum isn't 1, just normalize it.
So for an answer of Ans=[0.6, 0.0, 0.6] normalize it to [0.5, 0.0, 0.5]

Alan Coppola's method could work but your formulation of problem is wrong.

One one hand you say that weights should sum to 1 and on the other you say minimum weight is 40% here:

constraints=[sum_entries(w) == 1, w >= 0.40, w<=1.0]

hi, yes, i know the constraints are a contradiction and that is the whole point of me posting my question. I know the answer needs to be one weight at 0 and remaining 2 sum to 1. Will take a look at boolean and see if i can get an answer. If i target a portfolio based on 3 combinations, there will be one that has the lowest stdev and that will be the answer but how to code in python or which solver to use not sure. Will look at gurobi examples as well.

many thanks for the feedback so far

The constraints you've defined cannot be represented as a convex set. For CVXOPT or CVXPY, the feasible set must be convex. You can try MILP or something like that using Gurobi.

Quantopian supports:

http://docs.scipy.org/doc/scipy/reference/optimize.html

Have a look at:

http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html#scipy.optimize.minimize

I've been using:

Method SLSQP uses Sequential Least SQuares Programming to minimize a function of several variables with any combination of bounds, equality and inequality constraints. The method wraps the SLSQP Optimization subroutine originally implemented by Dieter Kraft [12]. Note that the wrapper handles infinite values in bounds by converting them into large floating values.

If you are interested, I could cook up an example. The bounds parameter can be used to limit the optimization search and the constraints can include the normalization, numpy.sum(x) = 1, where x is the asset allocation vector.

yes please, any example will be great and thank you for your time. For the three asset example

mu= [-0.01121522 -0.0029662 0.02089164]
covar = [[ 1.01555868e+00 2.82769766e-04 1.13099898e-03]
[ 2.82769766e-04 1.00342791e+00 -1.00797232e-02]
[ 1.13099898e-03 -1.00797232e-02 1.01961298e+00]]

minimize (1/2)x' covar x + mu' x where lower bound= 0.40 and upper =1 and sum(x) =1
will be cool to see how to set this up and see the weights.

thanks

Here's an example I had on-hand. It needs to be re-factored, but it should give you an idea of how things work. Note that SLSQP also accepts the Jacobian, so if an analytic form exists, you should probably include it, as well. It'll probably be several days before I can get back to this, but I think it'd be straightforward to tweak up the code I posted, with your objective function, etc. --Grant

610
Backtest from to with initial capital
Total Returns
--
Alpha
--
Beta
--
Sharpe
--
Sortino
--
Max Drawdown
--
Benchmark Returns
--
Volatility
--
 Returns 1 Month 3 Month 6 Month 12 Month
 Alpha 1 Month 3 Month 6 Month 12 Month
 Beta 1 Month 3 Month 6 Month 12 Month
 Sharpe 1 Month 3 Month 6 Month 12 Month
 Sortino 1 Month 3 Month 6 Month 12 Month
 Volatility 1 Month 3 Month 6 Month 12 Month
 Max Drawdown 1 Month 3 Month 6 Month 12 Month
# https://www.quantopian.com/posts/long-only-minimum-variance-portfolio-using-scipy-dot-optimize-dot-minimize

import numpy as np
from pytz import timezone
import scipy

def initialize(context):

context.stocks = [ sid(19662),  # XLY Consumer Discrectionary SPDR Fund
sid(19656),  # XLF Financial SPDR Fund
sid(19658),  # XLK Technology SPDR Fund
sid(19655),  # XLE Energy SPDR Fund
sid(19661),  # XLV Health Care SPRD Fund
sid(19657),  # XLI Industrial SPDR Fund
sid(19659),  # XLP Consumer Staples SPDR Fund
sid(19654),  # XLB Materials SPDR Fund
sid(19660),  # XLU Utilities SPRD Fund
sid(33652)]  # BND Vanguard Total Bond Market ETF

context.x0 = 1.0*np.ones_like(context.stocks)/len(context.stocks)

context.day_count = -1

def handle_data(context, data):

# Trade only once per day
loc_dt = get_datetime().astimezone(timezone('US/Eastern'))
if loc_dt.hour == 16 and loc_dt.minute == 0:
context.day_count += 1
pass
else:
return

if context.day_count % trading_freq != 0.0:
return

prices = history(21,'1d','price').as_matrix(context.stocks)
ret = np.diff(prices,axis=0) # daily returns

bnds = ((0,1),(0,1),(0,1),(0,1),(0,1),(0,1),(0,1),(0,1),(0,1),(0,1))
cons = ({'type': 'eq', 'fun': lambda x:  np.sum(x)-1.0})

res= scipy.optimize.minimize(variance, context.x0, args=ret,method='SLSQP',constraints=cons,bounds=bnds)

allocation = res.x
allocation[allocation<0]=0
denom = np.sum(allocation)
if denom != 0:
allocation = allocation/denom

context.x0 = allocation

record(stocks=np.sum(allocation[0:-1]))
record(bonds=allocation[-1])

for i,stock in enumerate(context.stocks):
order_target_percent(stock,allocation[i])

def variance(x,*args):

p = np.squeeze(np.asarray(args))
Acov = np.cov(p.T)

return np.dot(x,np.dot(Acov,x))
We have migrated this algorithm to work with a new version of the Quantopian API. The code is different than the original version, but the investment rationale of the algorithm has not changed. We've put everything you need to know here on one page.
There was a runtime error.

many thanks for your help....greatly appreciated.

You're welcome. Just to show that the bounds work, here's a version with limits on the bond portion of the portfolio.

Note that there are additional settings and you can see what the optimizer is doing, if you want:

http://docs.scipy.org/doc/scipy/reference/optimize.minimize-slsqp.html#optimize-minimize-slsqp

If you go this route, it is worth understanding the best settings. Also, you might try normalizing the returns to 1 in the example I posted, and then feeding them to the minimizer.

610
Backtest from to with initial capital
Total Returns
--
Alpha
--
Beta
--
Sharpe
--
Sortino
--
Max Drawdown
--
Benchmark Returns
--
Volatility
--
 Returns 1 Month 3 Month 6 Month 12 Month
 Alpha 1 Month 3 Month 6 Month 12 Month
 Beta 1 Month 3 Month 6 Month 12 Month
 Sharpe 1 Month 3 Month 6 Month 12 Month
 Sortino 1 Month 3 Month 6 Month 12 Month
 Volatility 1 Month 3 Month 6 Month 12 Month
 Max Drawdown 1 Month 3 Month 6 Month 12 Month
# https://www.quantopian.com/posts/long-only-minimum-variance-portfolio-using-scipy-dot-optimize-dot-minimize

import numpy as np
from pytz import timezone
import scipy

def initialize(context):

context.stocks = [ sid(19662),  # XLY Consumer Discrectionary SPDR Fund
sid(19656),  # XLF Financial SPDR Fund
sid(19658),  # XLK Technology SPDR Fund
sid(19655),  # XLE Energy SPDR Fund
sid(19661),  # XLV Health Care SPRD Fund
sid(19657),  # XLI Industrial SPDR Fund
sid(19659),  # XLP Consumer Staples SPDR Fund
sid(19654),  # XLB Materials SPDR Fund
sid(19660),  # XLU Utilities SPRD Fund
sid(33652)]  # BND Vanguard Total Bond Market ETF

context.x0 = 1.0*np.ones_like(context.stocks)/len(context.stocks)

context.day_count = -1

def handle_data(context, data):

# Trade only once per day
loc_dt = get_datetime().astimezone(timezone('US/Eastern'))
if loc_dt.hour == 16 and loc_dt.minute == 0:
context.day_count += 1
pass
else:
return

if context.day_count % trading_freq != 0.0:
return

prices = history(21,'1d','price').as_matrix(context.stocks)
ret = np.diff(prices,axis=0) # daily returns

bnds = ((0,1),(0,1),(0,1),(0,1),(0,1),(0,1),(0,1),(0,1),(0,1),(0.3,0.6))
cons = ({'type': 'eq', 'fun': lambda x:  np.sum(x)-1.0})

res= scipy.optimize.minimize(variance, context.x0, args=ret,method='SLSQP',constraints=cons,bounds=bnds)

allocation = res.x
allocation[allocation<0]=0
denom = np.sum(allocation)
if denom != 0:
allocation = allocation/denom

context.x0 = allocation

record(stocks=np.sum(allocation[0:-1]))
record(bonds=allocation[-1])

for i,stock in enumerate(context.stocks):
order_target_percent(stock,allocation[i])

def variance(x,*args):

p = np.squeeze(np.asarray(args))
Acov = np.cov(p.T)

return np.dot(x,np.dot(Acov,x))
We have migrated this algorithm to work with a new version of the Quantopian API. The code is different than the original version, but the investment rationale of the algorithm has not changed. We've put everything you need to know here on one page.
There was a runtime error.