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.

---

Maximizing Return

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

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\)

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

Formulation

The model maximizes the overall expected return for a prespecified maximum level of variance (risk). Mathematically, this results in a convex quadratically constrained optimization problem.

Decision Variables and Variable Bounds

The decision variables in the model are 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\).

Each position must be between 0 and 1; this prevents leverage and short-selling:

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

Constraints

The budget constraint ensures that all capital is invested:

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

The estimated risk must not exceed a prespecified maximal admissible level of variance \(\bar\sigma^2\):

\[x^\top \Sigma x \leq \bar\sigma^2\]

Objective Function

The objective is to maximize the expected return of the portfolio:

\[\max_x \mu^\top x\]

Using gurobipy, this can be expressed as follows:

V = 3.5  # maximal admissible variance (sigma^2)

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

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

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

# Limit on variance
risk_constr = m.addConstr(x @ Sigma.to_numpy() @ x <= V, name="Variance")

# Define objective function: Maximize expected return
m.setObjective(mu.to_numpy() @ x, gp.GRB.MAXIMIZE)

We now solve the optimization problem:

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 1 rows, 462 columns and 462 nonzeros
Model fingerprint: 0x1ba1fb96
Model has 1 quadratic constraint
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  QMatrix range    [3e-03, 1e+02]
  Objective range  [7e-02, 6e-01]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+00]
  QRHS range       [4e+00, 4e+00]
Presolve time: 0.05s
Presolved: 464 rows, 925 columns, 107878 nonzeros
Presolved model has 1 second-order cone constraint
Ordering time: 0.01s

Barrier statistics:
 AA' NZ     : 1.074e+05
 Factor NZ  : 1.079e+05 (roughly 1 MB of memory)
 Factor Ops : 3.341e+07 (less than 1 second per iteration)
 Threads    : 1

                  Objective                Residual
Iter       Primal          Dual         Primal    Dual     Compl     Time
   0   1.61202873e+01  1.74312373e-01  6.84e+01 6.67e-01  3.14e-02     0s
   1   2.06989937e+00  1.61037147e+00  7.63e+00 7.34e-07  4.32e-03     0s
   2   5.94842488e-01  8.44674441e-01  1.50e+00 1.41e-07  1.05e-03     0s
   3   2.41396673e-01  5.28985659e-01  1.65e-06 1.01e-08  2.07e-04     0s
   4   2.75348824e-01  4.18773720e-01  1.81e-12 4.33e-09  1.03e-04     0s
   5   3.17949843e-01  3.60990652e-01  6.15e-14 1.14e-09  3.10e-05     0s
   6   3.34976965e-01  3.43400301e-01  1.25e-13 9.17e-11  6.07e-06     0s
   7   3.37110375e-01  3.38160765e-01  6.08e-14 1.12e-11  7.57e-07     0s
   8   3.37545769e-01  3.37859496e-01  9.05e-14 1.41e-12  2.26e-07     0s
   9   3.37749001e-01  3.37776224e-01  1.95e-13 7.16e-14  1.96e-08     0s
  10   3.37772914e-01  3.37773790e-01  7.43e-12 2.14e-15  6.31e-10     0s

Barrier solved model in 10 iterations and 0.28 seconds (0.63 work units)
Optimal objective 3.37772914e-01

Display basic solution data:

print(f"Expected return:  {m.ObjVal:.6f}")
print(f"Variance:         {x.X @ Sigma @ x.X:.6f}")
print(f"Solution time:    {m.Runtime:.2f} seconds\n")

# Print investments (with non-negligible values, i.e., > 1e-5)
positions = pd.Series(name="Position", data=x.X, index=mu.index)
print(f"Number of trades: {positions[positions > 1e-5].count()}\n")
print(positions[positions > 1e-5])

Expected return:  0.337773
Variance:         3.499993
Solution time:    0.28 seconds

Number of trades: 25

TSLA    0.007032
KR      0.042558
PGR     0.126984
ORLY    0.048828
ODFL    0.032487
MNST    0.007642
KDP     0.081927
META    0.012825
UNH     0.017685
AVGO    0.048781
DXCM    0.012487
NFLX    0.015329
LLY     0.229288
DPZ     0.039872
MKTX    0.004647
WST     0.024237
TMUS    0.044854
NOC     0.048733
MOH     0.003070
MSFT    0.016673
WM      0.008693
TTWO    0.038648
ENPH    0.005489
NVDA    0.080938
AZO     0.000272
Name: Position, dtype: float64

Efficient Frontier

The efficient frontier reveals the balance between risk and return in investment portfolios. It shows the best-expected return level that can be achieved for a specified risk level. We compute this by solving the above optimization problem for a sample of admissible risk levels.

risks = np.linspace(1.37, 4, 20)
returns = np.zeros(risks.shape)
npos = np.zeros(risks.shape)

# hide Gurobi log output
m.params.OutputFlag = 0

# solve the model for each risk level
for i, risk_level in enumerate(risks):
    # set risk level: RHS of risk constraint
    risk_constr.QCRHS = risk_level**2
    m.optimize()
    # store data
    returns[i] = mu @ x.X
    npos[i] = len(x.X[x.X > 1e-5])

Next, we display the efficient frontier for this model: We plot the expected returns (on the \(y\)-axis) against the standard deviation \(\sqrt{x^\top\Sigma x}\) of the expected returns (on the \(x\)-axis). We also display the relationship between the risk and the number of positions in the optimal portfolio.

fig, axs = plt.subplots(1, 2, figsize=(10, 3))

# Axis 0: The efficient frontier
axs[0].scatter(x=risks, y=returns, marker="o", label="sample points", color="Red")
axs[0].plot(risks, returns, label="efficient frontier", color="Red")
axs[0].set_xlabel("Standard deviation")
axs[0].set_ylabel("Expected return")
axs[0].legend()
axs[0].grid()

# Axis 1: The number of open positions
axs[1].scatter(x=risks, y=npos, color="Red")
axs[1].plot(risks, npos, color="Red")
axs[1].set_xlabel("Standard deviation")
axs[1].set_ylabel("Number of positions")
axs[1].grid()

plt.show()

As expected, the number of open positions decreases as we allow more variance; the optimization will progressively invest in fewer high-risk but high-yield assets.