I am using python and the cvxopt library to calculate an efficient frontier, per the docs:
http://cvxopt.org/examples/book/portfolio.html
However, I cannot figure out how to add a constraint so that there is an upper bound on a particular asset's maximum allowed weight. Is that possible using cvxopt?
Here is my code so far that produces an efficient frontier with no constraints, except I believe b, which sets the max sum of weights to 1. I'm not sure what G, h, A, and mus do, and the docs don't really explain. Where does the 10**(5.0*t/N-1.0) in the formula for mus come from?
from math import sqrt
from cvxopt import matrix as cmatrix
from cvxopt.blas import dot as cdot
from cvxopt.solvers import qp, options
# Number of assets
n = 4
# Convariance matrix
S = cmatrix( [[ 4e-2, 6e-3, -4e-3, 0.0 ],
[ 6e-3, 1e-2, 0.0, 0.0 ],
[-4e-3, 0.0, 2.5e-3, 0.0 ],
[ 0.0, 0.0, 0.0, 0.0 ]] )
# Expected return
pbar = cmatrix([.12, .10, .07, .03])
# nxn matrix of 0s
G = cmatrix(0.0, (n,n))
# Convert G to negative identity matrix
G[::n+1] = -1.0
# nx1 matrix of 0s
h = cmatrix(0.0, (n,1))
# 1xn matrix of 1s
A = cmatrix(1.0, (1,n))
Aadd = cmatrix(0.0, (n,n))
# Convert Aadd to identity matrix
Aadd[::n+1] = 1.0
A = cmatrix(np.vstack((A,Aadd)))
# bounds
b = cmatrix([1.0,1.0,1.0,1.0,1.0])
N = 100
mus = [ 10**(5.0*t/N-1.0) for t in range(N) ]
options['show_progress'] = False
xs = [ qp(mu*S, -pbar, G, h, A, b)['x'] for mu in mus ]
returns = [ cdot(pbar,x) for x in xs ]
risks = [ sqrt(cdot(x, S*x)) for x in xs ]
sharpes = [x/y for x,y in zip(returns,risks)]
max_sharpe = max(sharpes)
max_index = sharpes.index(max_sharpe)
print xs[max_index]
Output:
File "<ipython-input-45-8e583df5adc5>", line 29, in <module>
xs = [ qp(mu*S, -pbar, G, h, A, b)['x'] for mu in mus ]
File "C:\Users\Anaconda\lib\site-packages\cvxopt\coneprog.py", line 4496, in qp
return coneqp(P, q, G, h, None, A, b, initvals)
File "C:\Users\Anaconda\lib\site-packages\cvxopt\coneprog.py", line 1986, in coneqp
raise ValueError("Rank(A) < p or Rank([P; G; A]) < n")
ValueError: Rank(A) < p or Rank([P; G; A]) < n