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.

[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 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:

[5]:

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:

[6]:

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 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.04s

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    : 2



                  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.90164131e-01  8.65923856e-01  1.48e+00 1.45e-07  1.05e-03     0s

   3   2.41151233e-01  5.13692754e-01  1.63e-06 9.67e-09  1.96e-04     0s

   4   2.81315045e-01  3.84260189e-01  1.80e-12 5.85e-10  7.42e-05     0s

   5   3.13602982e-01  3.66171811e-01  2.00e-13 1.19e-10  3.79e-05     0s

   6   3.34940007e-01  3.45092785e-01  1.12e-13 1.65e-11  7.31e-06     0s

   7   3.37712692e-01  3.38094058e-01  1.76e-13 7.59e-13  2.75e-07     0s

   8   3.37767816e-01  3.37775597e-01  3.38e-13 9.63e-15  5.61e-09     0s

   9   3.37772712e-01  3.37773661e-01  4.21e-12 2.97e-15  6.84e-10     0s



Barrier solved model in 9 iterations and 0.20 seconds (0.34 work units)

Optimal objective 3.37772712e-01

Display basic solution data:

[7]:

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.499989
Solution time:    0.20 seconds

Number of trades: 25

TSLA    0.006986
KR      0.042499
PGR     0.126965
ORLY    0.048806
ODFL    0.032537
MNST    0.007645
KDP     0.081846
META    0.012826
UNH     0.017542
AVGO    0.048707
DXCM    0.012490
NFLX    0.015361
LLY     0.229340
DPZ     0.039859
MKTX    0.004711
WST     0.024196
TMUS    0.044827
NOC     0.048677
MOH     0.003173
MSFT    0.016717
WM      0.008774
TTWO    0.038685
ENPH    0.005463
NVDA    0.080986
AZO     0.000366
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.

[8]:

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.

[9]:

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()
../_images/modeling_notebooks_basic_model_maxmu_15_0.png

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.