Optimization with turnover constraint

Question

I am optimizing using scipy.optimize using SLSQP. I am looking to minimize the variance with some upper bounds and lower bounds on each stock.

I am also looking to constraint the weight so that the turnover from the current portfolio is less than a certain limit. To do so, I have to include additional dummy variables (buy and sell for each stock).

Mathematically, I am trying to :

minimize $ w^T \Sigma w$

subject to $ \Sigma w = 1 $ and $ \Sigma \left|(w - w_o) \right| < c$

Since, absolute value will make the problem non-convex, I introduced two more set of variables, $w_b$ (buy weights) and $w_s$ (sell weights). And two additional constraints -

$ w_i = w_{o_i} + w_{b_i} + w_{s_i} $ for all stocks (This just means that the final weight is equal to starting weight + buy weight + sell weight)

$ w_{s_i} = w_{o_i} - w_{i}$ (Sell weight is forced equal to current weight minus final weight)

$ w_{b_i} = w_{i} - w_{o_i}$ (Buy weight is forced equal to final weight minus current weight)

Apart from this, there are other constraints like $0 < w_i < 0.3$ and $ max(0.3 - w_{o_i}, 0) > w_{b_i} => 0 $ and $ (-1 * w_{o_i}) < w_{s_i} <= 0 $

As of now, I haven't put the total turnover constraint.

Here's the toy code

import numpy as np
from scipy.optimize import minimize

def calculate_portfolio_var(w, V):
    if len(w) == V.shape[0]:
        return -1. * np.matmul(np.matmul(w.T, V), w)
    else:
        first_split = V.shape[0]
        w = w[: first_split]
        return -1. * np.matmul(np.matmul(w.T, V), w)

n_a = 10

par_wt = np.random.rand(n_a)
par_wt = par_wt / np.sum(par_wt)

cur_wt = np.random.rand(n_a)
cur_wt = cur_wt / np.sum(cur_wt)

ret = np.random.randn(n_a, 156)

lb = [0 for i in range(n_a * 3)]
ub = [0.3 for i in range(n_a * 3)]

lb[n_a * 2: ] =  -1 * cur_wt

ub[n_a: n_a * 2] = np.where(ub[:n_a] - cur_wt > 0, \
  ub[:n_a] - cur_wt, \
  np.repeat(0, n_a))

ub[n_a * 2: ] = np.where(lb[:n_a] - cur_wt > 0, \
  np.repeat(0, n_a), \
  lb[:n_a] - cur_wt > 0)

bnds = tuple(zip(lb, ub))

cov = np.cov(ret)

initial_guess = np.repeat(0.0, n_a * 3)
initial_guess[:n_a] = par_wt

cons = []

sum_cons = {'type' : 'eq', 
            'fun' : lambda x: np.sum(x) - 1}

cons.append(sum_cons)

for i in range(n_a):
    c1 = {'type' : 'eq', 
            'fun' : lambda x, i = i, cw = cur_wt, n = n_a : \
            (x[i] - x[n + i] + x[2*n + i] - cur_wt[i]) }

    c2 = {'type' : 'eq', 
            'fun' : lambda x, i = i, cw = cur_wt, n = n_a : \
            (x[n + i] - x[i] + cur_wt[i]) }

    c3 = {'type' : 'eq', 
            'fun' : lambda x, i = i, cw = cur_wt, n = n_a : \
            (x[2*n + i] - x[i] + cur_wt[i]) * 10}
    cons.append(c1)
    cons.append(c2)
    cons.append(c3)

res = minimize(calculate_portfolio_var, \
               args = cov, \
               bounds = bnds,\
               constraints = cons, \
               x0 = initial_guess, \
               method = 'SLSQP', \
               options = {'maxiter' : 10000, 'ftol' : 1e-15, \
                          'eps' : 1e-15, 'disp' : True})

print(res.status, np.round(res.fun, 3))

print(np.round(res.x[0:10], 3))

print(np.round(np.sum(res.x), 3))

print(np.round(np.sum(np.abs(res.x[:n_a] - cur_wt)), 4))

for i in range(10):    
    print(f'Current: {np.round(cur_wt[i], 4)}, Sell: {np.round(res.x[20+i], 4)}, Buy: {np.round(res.x[10+i], 4)}, Final: {np.round(res.x[i], 4)}')  

The output which I am getting is this. Obviously there is something wrong with this. Can someone please help me in formulating this in a way which scipy optimize understands?

More equality constraints than independent variables    (Exit mode 2)
            Current function value: -0.14473622934329722
            Iterations: 1
            Function evaluations: 32
            Gradient evaluations: 1
2 -0.145
[ 0.083  0.009  0.142  0.193  0.194  0.132  0.031  0.001  0.112  0.104]
1.0
0.8941
Current: 0.1718, Sell: 0.0, Buy: 0.0, Final: 0.0832
Current: 0.0511, Sell: 0.0, Buy: 0.0, Final: 0.0086
Current: 0.2313, Sell: 0.0, Buy: 0.0, Final: 0.1422
Current: 0.0589, Sell: 0.0, Buy: 0.0, Final: 0.1928
Current: 0.0003, Sell: 0.0, Buy: 0.0, Final: 0.1936
Current: 0.012, Sell: 0.0, Buy: 0.0, Final: 0.1318
Current: 0.0495, Sell: 0.0, Buy: 0.0, Final: 0.0306
Current: 0.0074, Sell: 0.0, Buy: 0.0, Final: 0.0013
Current: 0.213, Sell: 0.0, Buy: 0.0, Final: 0.1119
Current: 0.2047, Sell: 0.0, Buy: 0.0, Final: 0.104

Answer 1

You can probably also do this with scipy but there is a specific convex optimiser cvxopt that suits this problem.

from cvxopt import solvers, matrix

For your problem suppose you seek $N$ variables $w$, then let $x$ be a vector where the first $N$ elements are $w$ and the latter $N$ elements are slack variables $t$. i.e. $x=[w,t]^T$

For your problem $P = \begin{bmatrix} 2\Sigma & 0 \\ 0 & 0 \end{bmatrix}$ and $q = 0$:

N = 10 # size of problem
P = matrix(np.block([[2*Sigma, 0*Sigma], [0*Sigma, 0*Sigma]]))
q = matrix(np.zeros(2*N))

Subject to:

$$ \delta^T w = 1 $$ $$ \delta^T [|w-w_0|] \leq c $$

This can be reformulated in convex manner (https://scicomp.stackexchange.com/questions/20321/convex-optimization-problem-with-sum-of-absolute-value-constraints) as:

$$ [\delta^T, 0] x = 1 $$ $$ [0,\delta^T] x \leq c $$ $$ [-I_N, -I_N] x \leq -w_0 $$ $$ [I_N, -I_N] x \leq w_0 $$

A = matrix(np.block([np.ones(N), np.zeros(N)]))
b = matrix(1)
G = matrix(np.block([
       np.block([np.zeros(N), np.ones(N)]),
       np.block([-np.eye(N), -np.eye(N)]),
       np.block([np.eye(N), -np.eye(N)])
    ]))
h = matrix(np.block([[c], -w_0, w_0]))

soln = solvers.qp(P,q,G,h,A,b)

Note there may be a couple of bracket, or dimension bugs or sign issues, but this was a quick draft