# Volatility Target Optimization - Python

I'm testing a volatility target strategy in Python. This process involves solving the following optimization problem at each rebalance date

$$\min_w \left(w^T\Sigma w - \bar\sigma^2\right)^2$$ s.t. $$\mu^T w\ge \bar\mu$$ $$1^T w = 1$$

Which would be the best algorithm to implement this problem?

I'm assuming the barred variables are constants.

Clearly the minimum is obtained the closer the quadratic element gets to $\bar{\sigma}^2$ so the problem can be reformulated as:

$$\min_{w,t} t$$ s.t. $$w^T \Sigma w -t \leq \bar{\sigma}$$ $$w^T \Sigma w +t \geq \bar{\sigma}$$ $$\mu^T w \geq \bar{\mu}$$ $$1^Tw=1$$ $$t \geq 0$$

This is a Quadratically Constrained Linear Problem (QCLPs), which are a subclass of Quadratically Constrained Quadratic Programs (QCQPs), which are related to Second Order Cone Programs (SOCPs).

I don't know which algorithms will work best unfortunately but I would suggest Python's CVXOPT might have some documentation, or due to the nature of the problem I would hazard a guess that Sequential Quadratic Programming would be worth a try. Hopefully this answer can point you or someone else in the right direction.

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']