9.1 Markowitz portfolio optimization with chance constraints#
Preamble: Install Pyomo and a solver#
This cell selects and verifies a global SOLVER for the notebook. On Google Colab it installs Pyomo and Ipopt via the IDAES package and selects Ipopt. Elsewhere, it assumes Pyomo and a licensed Mosek installation are available and selects Mosek via the Pyomo SolverFactory.
The implementation below uses an explicit quadratic cone from pyomo.kernel. Mosek recognizes this cone directly; Ipopt handles its nonlinear representation. We initialize the cone variables away from the zero point and check solver termination before reading the solution.
import sys, os
if 'google.colab' in sys.modules:
%pip install idaes-pse --pre >/dev/null 2>/dev/null
!idaes get-extensions --to ./bin
os.environ['PATH'] += os.pathsep + os.path.abspath('bin')
solver = "ipopt"
else:
solver = "mosek_direct"
import pyomo.environ as pyo
SOLVER = pyo.SolverFactory(solver)
assert SOLVER.available(), f"Solver {solver} is not available."
import pyomo.kernel as pmo
import numpy as np
from scipy.stats import norm
Problem description#
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.
Suppose an initial capital \(C\) can be invested in \(n\) risky assets, each with an unknown return rate \(r_i\), \(i=1,\dots,n\). Let \(x_i\) be the amount invested in asset \(i\). When \(C=1\), these amounts are also fractions of the initial capital. 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}\) be 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 limit the probability that the return from the risky holdings falls below a threshold \(\alpha\). The risk-free contribution is included in the objective, but not in this chance constraint:
For positive variance, the risky return is normal with mean \(\mu^\top x\) and standard deviation \(\sqrt{x^\top\Sigma x}\). Hence,
Writing \(z_\beta=\Phi^{-1}(1-\beta)\), the chance constraint is equivalent to
The square root is essential: replacing the standard deviation by the variance changes the risk constraint. For \(0<\beta\leq 1/2\), \(z_\beta\) is nonnegative and the constraint is convex. The resulting portfolio optimization problem written as a SOCP is
We now implement this model for \(n=3\), \(\alpha=0.6\) and \(\beta=0.3\). For the positive definite covariance matrix used here, the Cholesky factorization gives \(\Sigma=B^\top B\), where B = np.linalg.cholesky(Sigma).T. Thus the chance constraint is the quadratic cone
We use the explicit cone interface introduced in Chapter 6, rather than passing a square-root expression to Mosek. The cone also enforces \(\mu^\top x-\alpha\geq0\); squaring the inequality without retaining this sign condition would not be equivalent.
This implementation requires positive definite \(\Sigma\). A singular covariance needs another factorization and separate treatment of zero-variance portfolios. Here \(\alpha>0\) forces every feasible portfolio to have positive variance.
# 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
# specify the initial capital, the risk-free return the number of risky assets,
# their expected returns, and their covariance matrix.
C = 1
R = 1.25
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]])
# To change covariance matrix, make sure you input a positive definite one.
# The easiest way to generate a random covariance matrix is first generating
# a full-rank m x m matrix A and then taking A^T A
# m = 3
# A = np.random.rand(m, m)
# Sigma = A.T @ A
def markowitz_chanceconstraints(alpha, beta, mu, Sigma, C=1, R=1.25):
if not 0 < beta <= 0.5:
raise ValueError("The convex formulation requires 0 < beta <= 0.5.")
n = len(mu)
if Sigma.shape != (n, n) or not np.allclose(Sigma, Sigma.T):
raise ValueError("Sigma must be a symmetric covariance matrix.")
# Positive definiteness is required by this Cholesky implementation.
B = np.linalg.cholesky(Sigma).T
z_beta = norm.ppf(1 - beta)
model = pmo.block()
model.x = pmo.variable_list(pmo.variable(lb=0, value=C / n) for _ in range(n))
model.xtilde = pmo.variable(lb=0, value=0)
model.risky_mean = pmo.expression(sum(mu[i] * model.x[i] for i in range(n)))
model.objective = pmo.objective(
model.risky_mean + R * model.xtilde, sense=pmo.maximize
)
model.total_assets = pmo.constraint(sum(model.x) + model.xtilde == C)
cone_coordinates = [
z_beta * sum(B[j, i] * model.x[i] for i in range(n)) for j in range(n)
]
model.chance_constraint = pmo.conic.quadratic.as_domain(
r=model.risky_mean - alpha, x=cone_coordinates
)
# Ipopt needs a nonzero starting point for the cone's nonlinear representation.
model.chance_constraint.r.value = max(1e-3, model.risky_mean() - alpha)
for variable, expression in zip(model.chance_constraint.x, cone_coordinates):
variable.value = pyo.value(expression)
return model
model = markowitz_chanceconstraints(alpha, beta, mu, Sigma, C=C, R=R)
result = SOLVER.solve(model)
pyo.assert_optimal_termination(result)
x = np.array([v.value for v in model.x])
risky_mean = float(mu @ x)
risky_std = float(np.sqrt(x @ Sigma @ x))
shortfall_probability = (
norm.cdf((alpha - risky_mean) / risky_std)
if risky_std > 0 else float(risky_mean <= alpha)
)
assert shortfall_probability <= beta + 1e-6
print(
f"Solver status: {result.solver.status}, Termination condition: {result.solver.termination_condition}"
)
print(f"Risk-free allocation: {model.xtilde.value:.6f}")
print("Risky allocations: " + ", ".join(f"{v:.6f}" for v in x))
print(f"Maximum expected return: {model.objective():.6f}")
print(f"Risky-return standard deviation: {risky_std:.6f}")
print(f"Shortfall probability: {shortfall_probability:.6f} (limit {beta:.6f})")
Solver status: ok, Termination condition: optimal
Risk-free allocation: 0.000000
Risky allocations: 0.666666, 0.116719, 0.216615
Maximum expected return: 1.259990
Risky-return standard deviation: 1.258560
Shortfall probability: 0.300000 (limit 0.300000)
For the stated data, the optimal allocation is approximately 0.6667, 0.1167 and 0.2166 in the risky assets, with no risk-free holding to the displayed precision. The maximum expected return is about 1.259990 per unit of initial capital. The shortfall probability is 0.30, so the risk limit is binding.
The earlier variance-based implementation gave an expected return of about 1.233250. Its actual shortfall probability was about 0.296342: it happened to be conservative for this instance, but replacing standard deviation by variance is not a valid general approximation.