Note

Download this Jupyter notebook and all data (unzip next to the ipynb file!). You will need a Gurobi license to run this notebook, please follow the license instructions.

Enforcing Diversification#

The standard mean-variance (Markowitz) portfolio selection model determines an optimal investment portfolio that balances risk and expected return. In this notebook, we minimize the variance (risk) of the portfolio, constraining the expected return to meet a prescribed minimum level. Please refer to the annotated list of references for more background information on portfolio optimization.

To this basic model, we add two types of diversification constraints:

  • Holdings must be diversified across a specified minimum number of assets.

  • Positions must fall between specified lower and upper bounds. This ensures that no single asset makes up more than a pre-specified percentage of the portfolio’s value.

[2]:

import gurobipy as gp
import pandas as pd
import numpy as np

Input Data#

The following input data is used within the model:

  • \(S\): set of stocks

  • \(\mu\): vector of expected returns

  • \(\Sigma\): PSD variance-covariance matrix

    • \(\sigma_{ij}\) covariance between returns of assets \(i\) and \(j\)

    • \(\sigma_{ii}\) variance of return of asset \(i\)

[4]:

# Import some example data set
Sigma = pd.read_pickle("sigma.pkl")
mu = pd.read_pickle("mu.pkl")

Formulation#

The model minimizes the variance of the portfolio given that the minimum level of expected return is attained and that some portfolio diversification requirements are met. Mathematically, this results in a convex quadratic mixed-integer optimization problem.

Model Parameters#

We use the following parameters:

  • \(\bar\mu\): required expected portfolio return

  • \(K\): minimal number of open positions in the portfolio

  • \(\ell > 0\): lower bound on position size

  • \(u \leq 1\): upper bound on position size

[5]:

# Values for the model parameters:
r = 0.25  # Required return
K = 35  # Minimal number of stocks
u = 0.15  # Maximal position size
l = 0.005  # Minimal position size

Decision Variables#

We need two sets of decision variables:

  1. The proportions of capital invested among the considered stocks. The corresponding vector of positions is denoted by \(x\) with its component \(x_i\) denoting the proportion of capital invested in stock \(i\).

  2. Binary variables \(b_i\) indicating whether or not asset \(i\) is held. If \(b_i\) is 0, the holding \(x_i\) is also 0; otherwise if \(b_i\) is 1, the investor holds asset \(i\) (that is, \(x_i \geq \ell\)).

Variable Bounds#

Each position \(x_i\) must be between 0 and \(u\); this prevents large positions as well as leverage and short-selling:

\begin{equation*} 0\leq x_i\leq u \; , \; i \in S\tag{1} \end{equation*}

The \(b_i\) must be binary:

\begin{equation*} b_i \in \{0,1\} \; , \; i \in S \end{equation*}

[6]:

%%capture
# Create an empty optimization model
m = gp.Model()

# Add variables: x[i] is the capital proportion invested in stock i; upper bound is u; see formula (1) above
x = m.addMVar(len(mu), lb=0, ub=u, name="x")

# Add binary variables: b[i]=1 if stock i is held; b[i]=0 otherwise
b = m.addMVar(len(mu), vtype=gp.GRB.BINARY, name="b")

Constraints#

The budget constraint ensures that all capital is invested:

\begin{equation*} \sum_{i \in S} x_i =1 \tag{2} \end{equation*}

The expected return of the portfolio must be at least \(\bar\mu\):

\begin{equation*} \mu^\top x \geq \bar\mu\tag{3} \end{equation*}

The variable bounds only enforce that each \(x_i\) is between \(0\) and \(u\). To enforce the minimal position size, we need the binary variables \(b\) and the following sets of discrete constraints:

Ensure that \(x_i = 0\) if \(b_i = 0\):

\begin{equation*} x_i \leq b_i \; , \; i \in S\tag{4} \end{equation*}

Note that \(x_i\) has an upper bound of \(u\leq 1\). Thus, if \(b_i = 1\), the above constraint is non-restrictive.

Ensure a minimal position size of \(\ell\) if asset \(i\) is traded:

\begin{equation*} x_i \geq \ell b_i \; , \; i \in S\tag{5} \end{equation*}

Hence \(b_i = 1\) implies \(x_i \geq \ell\). If \(b_i = 0\), this constraint is non-restrictive since \(x_i\) has a lower bound of 0.

Finally, there must be at least \(K\) positions in the portfolio:

\begin{equation*} \sum_{i \in S} b_i \geq K\tag{6} \end{equation*}

[7]:

%%capture
# Budget constraint: all investments sum up to 1; see formula (2) above
m.addConstr(x.sum() == 1, name="Budget_Constraint")

# Lower bound on expected return; see formula (3) above
m.addConstr(mu.to_numpy() @ x >= r, name="Minimal_Return")

# Force x to 0 if not traded; see formula (4) above
m.addConstr(x <= b, name="Indicator")

# Minimal position; see formula (5) above
m.addConstr(x >= l * b, name="Minimal_Position")

# Diversification constraint: at least B stocks must be held; see formula (6) above
m.addConstr(b.sum() >= K, name="Diversification")

Objective Function#

The objective is to minimize the risk of the portfolio, which is measured by its variance:

\[\min_x x^\top \Sigma x\]
[8]:

# Define objective function: Minimize risk
m.setObjective(x @ Sigma.to_numpy() @ x, gp.GRB.MINIMIZE)

We now solve the optimization problem:

[9]:

m.optimize()

Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (linux64 - "Ubuntu 22.04.4 LTS")



CPU model: Intel(R) Xeon(R) Platinum 8272CL CPU @ 2.60GHz, instruction set [SSE2|AVX|AVX2|AVX512]

Thread count: 2 physical cores, 2 logical processors, using up to 2 threads



WLS license 2443533 - registered to Gurobi GmbH

Optimize a model with 927 rows, 924 columns and 3234 nonzeros

Model fingerprint: 0x1c1570fe

Model has 106953 quadratic objective terms

Variable types: 462 continuous, 462 integer (462 binary)

Coefficient statistics:

  Matrix range     [5e-03, 1e+00]

  Objective range  [0e+00, 0e+00]

  QObjective range [6e-03, 2e+02]

  Bounds range     [1e-01, 1e+00]

  RHS range        [2e-01, 4e+01]

Found heuristic solution: objective 5.8417163

Presolve time: 0.04s

Presolved: 927 rows, 924 columns, 3233 nonzeros

Presolved model has 106953 quadratic objective terms

Variable types: 462 continuous, 462 integer (462 binary)



Root relaxation: objective 2.026834e+00, 1833 iterations, 0.03 seconds (0.05 work units)



    Nodes    |    Current Node    |     Objective Bounds      |     Work

 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time



     0     0    2.02683    0    8    5.84172    2.02683  65.3%     -    0s

H    0     0                       2.0423444    2.02683  0.76%     -    0s

     0     0    2.02695    0    8    2.04234    2.02695  0.75%     -    0s

H    0     0                       2.0270503    2.02695  0.01%     -    0s



Explored 1 nodes (1833 simplex iterations) in 0.27 seconds (0.12 work units)

Thread count was 2 (of 2 available processors)



Solution count 3: 2.02705 2.04234 5.84172



Optimal solution found (tolerance 1.00e-04)

Best objective 2.027050298073e+00, best bound 2.026945015707e+00, gap 0.0052%

Display basic solution data:

[10]:

print(f"Minimum Risk:     {m.ObjVal:.6f}")
print(f"Expected return:  {mu @ x.X:.6f}")
print(f"Solution time:    {m.Runtime:.2f} seconds\n")
print(f"Number of stocks: {sum(b.X)}\n")

# Print investments (with non-negligible value, i.e. >1e-5)
positions = pd.Series(name="Position", data=x.X, index=mu.index)
print(positions[positions > 1e-5])

Minimum Risk:     2.027050
Expected return:  0.250000
Solution time:    0.28 seconds

Number of stocks: 35.0

KR      0.029874
EG      0.005000
CTRA    0.005000
PGR     0.045229
CHRW    0.008562
CME     0.032737
ODFL    0.026182
BDX     0.024356
LIN     0.013054
MNST    0.005000
KDP     0.075960
GILD    0.016222
META    0.005006
CLX     0.055330
SJM     0.030111
PG      0.006777
LLY     0.102443
DPZ     0.052576
MKTX    0.019684
CPRT    0.005000
MRK     0.043846
ED      0.085151
WST     0.021548
TMUS    0.037665
NOC     0.016102
EA      0.005000
MSFT    0.005000
WM      0.044581
TTWO    0.036273
WMT     0.065286
TXN     0.005000
HRL     0.035113
XEL     0.005000
AZO     0.009048
CPB     0.021283
Name: Position, dtype: float64

Takeaways#

  • Cardinality constraints can be modeled using binary decision variables.

  • For diversification, it is essential to have non-zero lower bounds on the position size.