This blog post, as suggested by this answer, has a pretty complete implementation of a return maximize portfolio with target risk. The below is a simplified implementation I wrote for a previous project, for your reference. Note that the construction of matrices in cvxopt is not totally straightforward, and can take some getting used to when transitioning from np.
The referenced blog post uses mean historical return to create return forecasts, where my function below takes a return vector as input because I used different calculation method for forward-looking returns.
import cvxopt as cvx
def markowitz_opt(ret_vec, covar_mat, max_risk):
U,V = np.linalg.eig(covar_mat)
U[U<0] = 0
Usqrt = np.sqrt(U)
A = np.dot(np.diag(Usqrt), V.T)
# Calculating G and h matrix
G1temp = np.zeros((A.shape[0]+1, A.shape[1]))
G1temp[1:, :] = -A
h1temp = np.zeros((A.shape[0]+1, 1))
h1temp[0] = max_risk
ret_c = len(ret_vec)
for i in np.arange(ret_c):
ei = np.zeros((1, ret_c))
ei[0, i] = 1
if i == 0:
G2temp = [cvx.matrix(-ei)]
h2temp = [cvx.matrix(np.zeros((1,1)))]
else:
G2temp += [cvx.matrix(-ei)]
h2temp += [cvx.matrix(np.zeros((1,1)))]
# Construct list of matrices
Ftemp = np.ones((1, ret_c))
F = cvx.matrix(Ftemp)
g = cvx.matrix(np.ones((1,1)))
G = [cvx.matrix(G1temp)] + G2temp
H = [cvx.matrix(h1temp)] + h2temp
# Solce using QCQP
cvx.solvers.options['show_progress'] = False
sol = cvx.solvers.socp(
-cvx.matrix(ret_vec),
Gq=G, hq=H, A=F, b=g)
xsol = np.array(sol['x'])
return xsol, sol['status']