Markowitz portfolio optimization with chance constraints

Markowitz portfolio optimization with chance constraints#

# Install AMPL and solvers
%pip install -q amplpy pandas numpy scipy

SOLVER = "scip"

from amplpy import AMPL, ampl_notebook

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

?25l   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.0/5.7 MB ? eta -:--:--
   ━━━━╺━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.6/5.7 MB 17.4 MB/s eta 0:00:01
   ━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━━━━━━━━ 3.4/5.7 MB 49.4 MB/s eta 0:00:01
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸ 5.7/5.7 MB 64.2 MB/s eta 0:00:01
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 5.7/5.7 MB 47.1 MB/s eta 0:00:00
?25hUsing default Community Edition License for Colab. Get yours at: https://ampl.com/ce
Licensed to AMPL Community Edition License for the AMPL Model Colaboratory (https://ampl.com/colab).

from IPython.display import Markdown, HTML
import pandas as pd
import numpy as np

We consider here another variant of the Markowitz portfolio optimization problem, which we already encountered in the context of convex optimization here and in the context of conic optimization here.

Assuming there is an initial unit capital $C$ that needs to be invested in a selection of $n$ possible assets, each of them with a unknown return rate $r_i$, $i=1,\dots,n$. Let $x$ be the vector whose $i$-th component $x_i$ describes the fraction of the capital invested in asset $i$. The return rate vector $r$ can be modelled by a multivariate Gaussian distribution with mean $\mu$ and covariance $\Sigma$. Assume there is also a risk-free asset with guaranteed return rate $R$ and let $\tilde{x}$ the amount invested in that asset. We want to determine the portfolio that maximizes the expected return $\mathbb{E} ( R \tilde{x} + r^\top x )$, which in view of our assumptions rewrites as $ \mathbb{E} ( R \tilde{x} + r^\top x ) = R \tilde{x} + \mu^\top x$.

Additionally, we includ a loss risk chance constraint of the form

\[ \mathbb{P} ( r^\top x \leq \alpha) \leq \beta. \]

Thanks to the assumption that $r$ is multivariate Gaussian, this chance constraint can be equivalently rewritten as

\[ \mu^\top x \geq \Phi^{-1}(1-\beta) \| \Sigma^{1/2} r \|_2 + \alpha, \]

which the theory guarantees to be a convex constraint if $\beta \leq 1/2$. The resulting portfolio optimization problem written as a SOCP is

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

We now implement an AMPL model and solve it for $n=3$, $\alpha = 0.6$, $\beta =0.3$, the given vector $\mu$ and semi-definite positive covariance matrix $\Sigma$.

# we import the inverse CDF or quantile function for the standard normal norm.ppf() from scipy.stats
from scipy.stats import norm

# We set our risk threshold and risk levels (sometimes you may get an infeasible problem if the chance
# constraint becomes too tight!)
alpha = 0.6
beta = 0.3

# We specify the initial capital, the risk-free return the number of risky assets, their expected returns, and their covariance matrix.
C = 1
R = 1.05
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]])

# Check how dramatically the optimal solution changes if we assume i.i.d. deviations for the returns.
# Sigma = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])

# If you want to change covariance matrix, make sure 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

%%writefile markowitz_chanceconstraints.mod

set N;

param Sigma{N, N};
param mu{N};
param C;
param R;
param gamma;

param phi_val;
param alpha;

var x_tilde >= 0;
var x{N} >= 0;

maximize expected_return:
    sum{i in N} mu[i] * x[i];

subject to bounded_variance:
    phi_val * (sum{i in N, j in N} x[i] * Sigma[i, j] * x[j]) <= sum{i in N} mu[i] * x[i] - alpha;

subject to total_assets:
    sum{i in N} x[i] + x_tilde = C;

Writing markowitz_chanceconstraints.mod

def markowitz_chanceconstraints(alpha, beta, mu, Sigma):
    m = AMPL()
    m.read("markowitz_chanceconstraints.mod")

    Sigma_df = pd.DataFrame(Sigma, index=range(n), columns=range(n))
    phi_val = norm.ppf(1 - beta)

    m.set["N"] = range(n)

    m.param["phi_val"] = phi_val
    m.param["alpha"] = alpha
    m.param["C"] = C
    m.param["R"] = R
    m.param["mu"] = mu
    m.param["Sigma"] = Sigma_df

    m.solve(solver=SOLVER)
    return m.solve_result, m


solve_result, model = markowitz_chanceconstraints(alpha, beta, mu, Sigma)
assert solve_result == "solved", solve_result

x_tilde = model.var["x_tilde"].value()
x = model.var["x"].to_dict()

display(Markdown(f"**Solver result:** *{solve_result}*"))
display(
    Markdown(
        f"**Solution:** $\\tilde x = {x_tilde:.3f}$, $x_1 = {x[0]:.3f}$,  $x_2 = {x[1]:.3f}$,  $x_3 = {x[2]:.3f}$"
    )
)
display(
    Markdown(
        f"**Maximizes objective value to:** ${model.obj['expected_return'].value():.2f}$"
    )
)

SCIP 9.0.1: WARNING: No dual information available when presolving was performed.
SCIP 9.0.1: optimal solution; objective 1.232266138
685 simplex iterations

Solver result: solved

Solution: $\tilde x = 0.000$, $x_1 = 0.666$, $x_2 = 0.256$, $x_3 = 0.078$

Maximizes objective value to: $1.23$