Markowitz portfolio optimization revisited

# install dependencies and select solver
%pip install -q amplpy numpy pandas matplotlib

SOLVER_CONIC = "mosek"  # ipopt, mosek, gurobi, knitro

from amplpy import AMPL, ampl_notebook

ampl = ampl_notebook(
    modules=["coin", "mosek"],  # modules to install
    license_uuid="default",  # license to use
)  # instantiate AMPL object and register notebook magics

from IPython.display import Markdown, HTML
import numpy as np
import matplotlib.pyplot as plt

Problem description and model formulation

Consider again the Markowitz portfolio optimization we presented earlier in Chapter 5. Recall that the matrix \(\Sigma\) describes the covariance among the uncertain return rates \(r_i\), \(i=1,\dots, n\). Since \(\Sigma\) is positive semidefinite by definition, it allows for a Cholesky factorization, namely \(\Sigma = B B^\top\). We can then rewrite the quadratic constraint as \(\|B^\top x \|_2 \leq \gamma\) and thus as \((\gamma, B^\top x) \in \mathcal{L}^{n+1}\) using the Lorentz cone. In this way, we realize that the original portfolio problem we formulated earlier is in fact a conic quadratic optimization problem, which can thus be solved faster and more reliably. The optimal solution of that problem was the one with the maximum expected return while allowing for a specific level \(\gamma\) of risk.

However, an investor could aim for a different trade-off between return and risk and formulate a slightly different optimization problem, namely

\[\begin{split} \begin{align*} \max \quad & R \tilde{x} + \mu^\top x - \alpha x^\top \Sigma x \\ \text{s.t.}\quad & \sum_{i=1}^n x_i + \tilde{x} = C \\ & \tilde x \geq 0\\ & x_i \geq 0 & \forall \, i=1,\dots,n. \end{align*} \end{split}\]

where \(\alpha \geq 0\) is a risk tolerance parameter that describes the relative importance of return vs. risk for the investor.

The risk, quantified by the variance of the investment return \(x^\top \Sigma x = x^\top B^\top B x\), appears now in the objective function as a penalty term. Note that even in this new formulation we have a conic problem since we can rewrite it as

\[\begin{split} \begin{align*} \max \quad & R \tilde{x} + \mu^\top x - \alpha s \\ \text{s.t.}\quad & \sum_{i=1}^n x_i + \tilde{x} = C \\ & \| B^\top x\|^2_2 \leq s \\ & \tilde x \geq 0 \\ & s \geq 0\\ & x_i \geq 0 & \forall \, i=1,\dots,n. \end{align*} \end{split}\]

Solving for all values of \(\alpha \geq 0\), one can obtain the so-called efficient frontier.

%%writefile markowitz_portfolio_revisited.mod

# Specify the initial capital, the risk tolerance, and the guaranteed return rate. 
param C;
param alpha;
param R;

# Specify the number of assets, their expected return, and their covariance matrix.
param n;
param mu{1..n};
param Sigma{1..n, 1..n};

var xtilde >= 0;
var x{1..n} >= 0;
var s >= 0;
    
maximize Objective:
    sum {i in 1..n} mu[i]*x[i] + R*xtilde - alpha*s;
    
s.t. BoundedVariance:
    sum {i in 1..n, j in 1..n} x[i]*Sigma[i,j]*x[j] <= s;
    
s.t. TotalAssets:
    sum {i in 1..n} x[i] + xtilde == C;
    

Overwriting markowitz_portfolio_revisited.mod

# Specify the initial capital, the risk tolerance, and the guaranteed return rate.
C = 1
alpha = 0.1
R = 1.05

# Specify the number of assets, their expected return, and their covariance matrix.
n = 3
mu = np.array([1.25, 1.15, 1.35])
Sigma = np.array([[1.5, 0.5, 2], [0.5, 2, 0], [2, 0, 5]])

# If you want to change the covariance matrix Sigma, ensure you input a semi-definite positive one.
# The easiest way to generate a random covariance matrix is first generating a random m x m matrix A
# and then taking the matrix A^T A (which is always semi-definite positive)
# m = 3
# A = np.random.rand(m, m)
# Sigma = A.T @ A
#
# Moreover, in practive such a matrix A, called factor, can be low-rank,
# see https://docs.mosek.com/modeling-cookbook/qcqo.html#example-factor-model.
# This would provide better numerical properties for the proper conic formulation
#        y=Ax, |y|^2 <= s,
# corresponding to the mathematical formulation above.


def markowitz_revisited(alpha, mu, Sigma):
    ampl = AMPL()
    ampl.read("markowitz_portfolio_revisited.mod")

    ampl.param["C"] = C
    ampl.param["alpha"] = alpha
    ampl.param["R"] = R

    ampl.param["n"] = n
    ampl.param["mu"] = mu
    ampl.param["Sigma"] = {
        (i + 1, j + 1): Sigma[i][j] for i in range(n) for j in range(n)
    }

    ampl.solve(solver=SOLVER_CONIC)
    assert ampl.solve_result == "solved", ampl.solve_result

    return ampl


ampl = markowitz_revisited(alpha, mu, Sigma)
xtilde = ampl.get_variable("xtilde").value()
x = ampl.get_variable("x").to_dict()

display(Markdown(f"**Solver status:** *{ampl.solve_result}*"))
display(
    Markdown(
        f"**Solution:** $\\tilde x = {xtilde:.3f}$"
        f"              $x_1 = {x[1]:.3f}$,  $x_2 = {x[2]:.3f}$,  $x_3 = {x[3]:.3f}$"
    )
)
display(
    Markdown(
        f"**Maximizes objective value to:**"
        f" ${ampl.get_objective('Objective').value():.2f}$"
    )
)

MOSEK 10.0.43: MOSEK 10.0.43: optimal; objective 1.12065217
0 simplex iterations
8 barrier iterations

Solver status: solved

Solution: $\tilde x = 0.283$ $x_1 = 0.478$, $x_2 = 0.130$, $x_3 = 0.109$

Maximizes objective value to: $1.12$

alpha_values = [
    0.005,
    0.01,
    0.02,
    0.03,
    0.04,
    0.05,
    0.06,
    0.07,
    0.08,
    0.09,
    0.1,
    0.11,
    0.12,
    0.13,
    0.14,
    0.15,
    0.16,
    0.17,
    0.18,
    0.19,
    0.20,
    0.25,
    0.5,
]
objective = []

for alpha in alpha_values:
    ampl = markowitz_revisited(alpha, mu, Sigma)
    objective.append(round(ampl.get_objective("Objective").value(), 3))

plt.plot(alpha_values, objective)
plt.xlabel(r"Risk tolerance $\alpha$")
plt.ylabel("Optimal objective value")
plt.show()

MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.324999981
0 simplex iterations
8 barrier iterations
MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.299999995
0 simplex iterations
8 barrier iterations
MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.251999955
0 simplex iterations
9 barrier iterations
MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.221333307
0 simplex iterations
8 barrier iterations
MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.199047619
0 simplex iterations
10 barrier iterations
MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.180952352
0 simplex iterations
7 barrier iterations
MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.165238086
0 simplex iterations
8 barrier iterations
MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.150884338
0 simplex iterations
9 barrier iterations
MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.138315214
0 simplex iterations
7 barrier iterations
MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.128502403
0 simplex iterations
8 barrier iterations
MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.12065217
0 simplex iterations
8 barrier iterations
MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.114229246
0 simplex iterations
8 barrier iterations
MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.108876806
0 simplex iterations
8 barrier iterations
MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.104347815
0 simplex iterations
8 barrier iterations
MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.100465835
0 simplex iterations
8 barrier iterations
MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.097101441
0 simplex iterations
8 barrier iterations
MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.094157586
0 simplex iterations
8 barrier iterations
MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.091560078
0 simplex iterations
8 barrier iterations
MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.089251187
0 simplex iterations
8 barrier iterations
MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.087185341
0 simplex iterations
8 barrier iterations
MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.085326081
0 simplex iterations
8 barrier iterations
MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.078260869
0 simplex iterations
8 barrier iterations
MOSEK 10.0.43:MOSEK 10.0.43: optimal; objective 1.064130422
0 simplex iterations
9 barrier iterations