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 Utility#

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 utility, which is defined as a weighted linear function of return and risk that can be adjusted by varying the risk-aversion coefficient \(\gamma\). 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 investor’s expected utility, which is defined as a linear combination of the expected return and risk, accounting for the investor’s level of risk aversion via the risk-aversion coefficient \(\gamma\). Mathematically, this results in a convex quadratic 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\]

Objective Function#

The objective is to maximize expected utility:

\[\max_x \mu^\top x - \tfrac{\gamma}{2} x^\top \Sigma x\]

Using gurobipy, this can be expressed as follows:

[5]:

%%capture
gamma = 0.2  # risk-aversion coefficient

# 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")

# Define objective function: Maximize expected utility
m.setObjective(
    mu.to_numpy() @ x - gamma / 2 * (x @ Sigma.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: 0xee0d2d3c

Model has 106953 quadratic objective terms

Coefficient statistics:

  Matrix range     [1e+00, 1e+00]

  Objective range  [7e-02, 6e-01]

  QObjective range [6e-04, 2e+01]

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

  RHS range        [1e+00, 1e+00]

Presolve time: 0.03s

Presolved: 1 rows, 462 columns, 462 nonzeros

Presolved model has 106953 quadratic objective terms

Ordering time: 0.01s



Barrier statistics:

 Free vars  : 461

 AA' NZ     : 1.065e+05

 Factor NZ  : 1.070e+05 (roughly 1 MB of memory)

 Factor Ops : 3.298e+07 (less than 1 second per iteration)

 Threads    : 2



                  Objective                Residual

Iter       Primal          Dual         Primal    Dual     Compl     Time

   0  -2.00616666e+05  3.17727204e+05  2.52e+05 6.60e-01  2.50e+05     0s

   1  -1.03370221e+05  1.09628139e+05  1.26e+04 1.69e-05  1.32e+04     0s

   2  -6.68666276e+02  1.19323021e+03  1.10e+02 1.47e-07  1.20e+02     0s

   3  -5.68494474e-01  4.96535897e+02  1.10e-04 1.47e-13  5.38e-01     0s

   4  -5.66713960e-01  1.32387974e+00  2.93e-07 1.46e-13  2.05e-03     0s

   5  -1.97476330e-01  3.51421951e-01  1.19e-08 5.77e-15  5.94e-04     0s

   6  -3.92447116e-02  1.81471841e-01  1.32e-14 5.00e-16  2.39e-04     0s

   7   1.65129388e-02  1.17129665e-01  5.55e-16 3.89e-16  1.09e-04     0s

   8   3.96401803e-02  5.48244228e-02  5.55e-16 2.22e-16  1.64e-05     0s

   9   4.76945608e-02  5.03333839e-02  1.66e-15 4.44e-16  2.86e-06     0s

  10   4.98128643e-02  4.98632128e-02  2.36e-15 3.33e-16  5.45e-08     0s

  11   4.98587716e-02  4.98591706e-02  2.52e-14 2.78e-16  4.32e-10     0s

  12   4.98590054e-02  4.98590288e-02  2.20e-11 2.78e-16  2.54e-11     0s

  13   4.98590147e-02  4.98590163e-02  3.98e-11 2.78e-16  1.76e-12     0s



Barrier solved model in 13 iterations and 0.27 seconds (0.45 work units)

Optimal objective 4.98590147e-02

Display basic solution data:

[7]:

print(f"Expected utility: {m.ObjVal:.6f}")
print(f"Minimum risk:     {x.X @ Sigma @ x.X:.6f}")
print(f"Expected return:  {mu @ 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 utility: 0.049859
Minimum risk:     2.158453
Expected return:  0.265704
Solution time:    0.28 seconds

Number of trades: 31

KR      0.034830
PGR     0.065163
CME     0.028443
ODFL    0.036042
BDX     0.017070
MNST    0.007179
BR      0.000587
KDP     0.082628
GILD    0.007109
META    0.009799
CLX     0.046160
SJM     0.016283
PG      0.000170
LLY     0.128341
DPZ     0.054559
MKTX    0.019686
CPRT    0.003111
MRK     0.030565
ED      0.075465
WST     0.023673
TMUS    0.045910
NOC     0.024009
EA      0.001504
MSFT    0.012610
WM      0.050315
TTWO    0.041693
WMT     0.058442
NVDA    0.009022
HRL     0.031237
AZO     0.019367
CPB     0.019028
Name: Position, dtype: float64

Efficient Frontier#

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

[8]:

# Hand picked for approximately equidistant return/risk points
# fmt: off
gamma_vals = np.array([
    2.00e-03,3.10e-03,9.67e-03,1.10e-02,1.17e-02,1.26e-02,1.35e-02,1.47e-02,1.62e-02,1.81e-02,
    2.05e-02,2.36e-02,2.69e-02,3.12e-02,3.69e-02,4.39e-02,5.52e-02,7.36e-02,1.10e-01,2.21e-01,])
# fmt: on

returns = np.zeros(gamma_vals.shape)
risks = np.zeros(gamma_vals.shape)
npos = np.zeros(gamma_vals.shape)

utility = []
# prevent Gurobi log output
m.params.OutputFlag = 0

# solve the model for each risk level
for i, gamma in enumerate(gamma_vals):
    # Replace utility objective function according to this choice for gamma
    m.setObjective(
        mu.to_numpy() @ x - gamma / 2.0 * (x @ Sigma.to_numpy() @ x), gp.GRB.MAXIMIZE
    )
    m.optimize()
    utility.append(m.objval)
    # store data
    returns[i] = mu @ x.X
    risks[i] = np.sqrt(x.X @ Sigma @ 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_maxutility_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.