How to calculate optimal portfolio using sector constraints in python

Question

I'm looking into CVXPY at the moment.

Main goal would be to be able to calculate the optimal portfolio, which in my opinion would mean that we need to maximise

(expected return - risk free) / volatility

To make it simpler I would like to drop risk free out of the equation (it's anyway near zero at the moment).

My problem is: I would like to constrain the weights in the portfolio based on maximum allocation that can be allocated to a certain sector.

Sector   Stocks
A        1
         2
         3
B        4
         5
         6

how can i achieve, that the minimum and maximum of an allocation to every sector lies in between some bands?

I've already seen this post: industry level constraints

But i really would like to implement it with cvxpy as the answer suggested. Unfortunately I have no idea how to implement this?

I already got the covarmatrix, return estimates and bands on hand but simply can't find a way how to implement those constraints. Any help would be highly appreciated!

EDIT: I tried my very best and came up with this script: thanks to the comment from @Attack68 in combination with this paper: Sharpe Quadratic Optimization

import cvxpy as cp
import numpy as np 

np.random.seed(101)

## NUMBER OF ASSETS
n_assets = 4

## NUMBER OF OBSERVATIONS
n_obs = 1000

## GENERATE RANDOM RETURNS
return_vec = np.random.randn(n_assets, n_obs)

## SET UP PROBLEM
C = np.asmatrix(np.cov(return_vec))          # Covar Matrix
mu = np.asmatrix(np.mean(return_vec,axis=1)) # return estimat
mu0 = -0.0075                                # risk free rate

y = cp.Variable(n_assets)                    # "weights"

A = np.asmatrix([[0.6,0.6,0,0],[-1.2,-1.2,0,0],[0,0,0.2,0.2],[0,0,-1.2,-1.2]])
bounds = np.asmatrix([0.4,0.2,0.8,0.2])

# HOW CAN I SUBSTRACT THE BOUNDS FROM ONLY THE NON-NULL VALUES? (AS FAR AS THIS WILL BE NEEDED?)
A_mod = A #- bounds.T


## CREATE CONSTRAINTS
constraints = [(mu-mu0)@y==1,
               y >= 0,
               y@A_mod >= 0]

## FORM OBJECTIVE
obj = cp.Minimize(cp.quad_form(y,C))

## FORM AND SOLVE PROBLEM
prob = cp.Problem(obj, constraints)
prob.solve()
w = y.value/sum(y.value)
w

array([0.35386785, 0.17693393, 0.4034641 , 0.06573412])

Unfortunately this solution doesn't meet the constraints. What am I doing wrong?

Answer 1

This is very simply done. It involves ensuring the constraints are presented as part of the matrix standard form.

You will typically have the constraint that all assets sum to one, i.e. the matrix-vector equation:

$$ \delta^T x = 1 $$

If you want to create an inequality constraint for assets in a sector just isolate them:

$$ \begin{bmatrix} 1 & 1 & 1 & 0 & 0 & 0 \\ -1 & -1 & -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 1 & 1 \\ 0 & 0 & 0 & -1 & -1 & -1 \\ \end{bmatrix} x \leq \begin{bmatrix} a_1 \\ a_2 \\ b_1 \\ b_2 \end{bmatrix}$$

Then add(stack) these matrix-vector rows to any existing inequality constraints you have.

Answer 2

For those experiencing a similar problem, here is the solution that worked for me:

## OPTIMIZE PORTFOLIO WEIGHTS UNDER THE OBJECTIVE OF MAXIMIZING THE SHARPE RATIO
## WHILE CONSTRAINING THE WEIGHTS TO SECTOR BOUNDS
## PAPER: # ACCORDING TO: http://people.stat.sc.edu/sshen/events/backtesting/reference/maximizing%20the%20sharpe%20ratio.pdf

np.random.seed(101)

## NUMBER OF ASSETS
n_assets = 4

## NUMBER OF OBSERVATIONS
n_obs = 1000

## GENERATE RANDOM RETURNS
return_vec = np.random.randn(n_assets, n_obs)

## SET UP PROBLEM
C = np.asmatrix(np.cov(return_vec))
mu = np.asmatrix(np.mean(return_vec,axis=1))
mu0 = np.asmatrix(-0.0075)

## INITIATE WEIGHT VARIABLE
y = cp.Variable(n_assets)

# DEFINE CONSTRAINTS AND MODIFY FOR QUADRATIC PROBLEM
A = np.asmatrix([[1,1,0,0],[-1,-1,0,0],[0,0,1,1],[0,0,-1,-1]])
bounds = np.asmatrix([0.2,-0.5,0.2,-0.8])
A_mod = A - bounds.T

## CREATE CONSTRAINTS
constraints = [(mu-mu0)@y==1,
               y >= 0,
               [email protected] >= 0]

## FORM OBJECTIVE
obj = cp.Minimize(cp.quad_form(y,C))

## FORM AND SOLVE PROBLEM
prob = cp.Problem(obj, constraints)
prob.solve()

## TRANSFORM FINAL WEIGHTS
w = y.value/sum(y.value)