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.


Cardinality Constraints

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 cardinality constraints to

  • limit the number of open positions and

  • limit the minimum size of each position.

These limits are often used to reduce transaction costs and simplify the rebalancing and management of the portfolio.

[2]:
import gurobipy as gp
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

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, that the number of positions held does not exceed a specified number, and that the size of each position does not fall below a certain level.

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\): maximal number of stocks in the portfolio

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

[5]:
# Values for the model parameters:
r = 0.25  # Required return
K = 15  # Maximal number of stocks
l = 0.00001  # 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 must be between 0 and 1; this prevents leverage and short-selling:

\[0\leq x_i\leq 1 \; , \; i \in S\]

The \(b_i\) must be binary:

\[b_i \in \{0,1\} \; , \; i \in S\]
[6]:
%%capture
# Create an empty optimization model
m = gp.Model()

# Add variables: x[i] denotes the proportion invested in stock i
x = m.addMVar(len(mu), lb=0, ub=1, name="x")

# Add variables: b[i]=1 if stock i is held, and 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:

\[\sum_{i \in S} x_i =1\]

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

\[\mu^\top x \geq \bar\mu\]

The variable bounds only enforce that each \(x_i\) is between \(0\) and \(1\). To enforce the minimal position size, we use 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{1} \end{equation*}

Note that \(x_i\) has an upper bound of 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{2} \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 most \(K\) positions in the portfolio:

\begin{equation*} \sum_{i \in S} b_i \leq K \end{equation*}

[7]:
# Budget constraint: all investments sum up to 1
m.addConstr(x.sum() == 1, "Budget_Constraint")

# Lower bound on expected return
m.addConstr(mu.to_numpy() @ x >= r, "Minimal_Return")

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

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

# Cardinality constraint: at most K positions
cardinality_constr = m.addConstr(b.sum() <= K, "Cardinality")

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: AMD EPYC 7763 64-Core Processor, instruction set [SSE2|AVX|AVX2]
Thread count: 1 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: 0x73f3de97
Model has 106953 quadratic objective terms
Variable types: 462 continuous, 462 integer (462 binary)
Coefficient statistics:
  Matrix range     [1e-05, 1e+00]
  Objective range  [0e+00, 0e+00]
  QObjective range [6e-03, 2e+02]
  Bounds range     [1e+00, 1e+00]
  RHS range        [2e-01, 2e+01]
Found heuristic solution: objective 21.1318982
Presolve time: 0.05s
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.026597e+00, 154 iterations, 0.01 seconds (0.01 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0    2.02660    0   37   21.13190    2.02660  90.4%     -    0s
H    0     0                       2.5684442    2.02660  21.1%     -    0s
     0     0    2.02660    0   37    2.56844    2.02660  21.1%     -    0s
H    0     0                       2.0761887    2.02660  2.39%     -    0s
     0     0    2.02660    0   37    2.07619    2.02660  2.39%     -    0s
     0     0    2.02660    0   37    2.07619    2.02660  2.39%     -    0s
     0     0    2.02660    0   37    2.07619    2.02660  2.39%     -    0s
     0     0    2.02660    0   36    2.07619    2.02660  2.39%     -    0s
H    0     0                       2.0717477    2.02660  2.18%     -    0s
     0     2    2.02660    0   36    2.07175    2.02660  2.18%     -    0s
*  130    51              18       2.0661170    2.03162  1.67%   7.5    0s
  4590   880    2.06573   28   21    2.06612    2.04366  1.09%   6.3    5s
 10226  1247    2.06178   29   23    2.06612    2.05328  0.62%   5.9   10s
 15980   977    2.06202   29   24    2.06612    2.05902  0.34%   5.9   15s

Cutting planes:
  Implied bound: 5
  MIR: 12
  Flow cover: 1

Explored 20390 nodes (117668 simplex iterations) in 18.15 seconds (27.73 work units)
Thread count was 2 (of 2 available processors)

Solution count 5: 2.06612 2.07175 2.07619 ... 21.1319

Optimal solution found (tolerance 1.00e-04)
Best objective 2.066116973868e+00, best bound 2.065910672679e+00, gap 0.0100%

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 trades: {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.066117
Expected return:  0.250000
Solution time:    18.15 seconds

Number of trades: 15.0

KR      0.040359
PGR     0.051136
CME     0.049374
ODFL    0.044432
KDP     0.087797
CLX     0.073365
SJM     0.060471
LLY     0.104751
DPZ     0.058761
MRK     0.073922
ED      0.110132
TMUS    0.047234
WM      0.073017
TTWO    0.047324
WMT     0.077925
Name: Position, dtype: float64

Comparison with the unconstrained portfolio

We can also compute the portfolio without the cardinality constraint and compare the resulting portfolios.

[11]:
# remove cardinality constraint
m.remove(cardinality_constr)
m.params.OutputFlag = 0
m.optimize()

# retrieve and display solution data
positions_unconstr = pd.Series(name="Position", data=x.X, index=mu.index)
mask = (positions > 1e-5) | (positions_unconstr > 1e-5)
df = pd.DataFrame(
    index=positions[mask].index,
    data={
        "at most 15 assets": positions,
        "unconstrained": positions_unconstr,
    },
).sort_values(by=["at most 15 assets", "unconstrained"], ascending=True)

axs = df.plot.barh(color=["#0b1a3c", "#dd2113"])
axs.set_xlabel("Fraction of investment sum")
plt.title("Minimum Variance portfolios with and without cardinality constraint")
plt.show()
../_images/modeling_notebooks_concentration_19_0.png

Takeaways

  • Cardinality constraints can be modeled using binary decision variables.

  • This already leads to non-trivial combinatorial optimization problems, which can be solved using Gurobi.