Solving a system of two equations with non-convex matrix multiplication for MV optimization

Question

Scenario: I am trying to do a variation of the MV optimization for a portfolio. In this instance, I already have a vector of mean returns ($\mu$), a vector of ones, a covariance matrix ($\Sigma$), and $\phi^{-1}$ which is the inverse of the standard normal cumulative distribution function (function of probability alpha).

Problem: From the two equations stated below, I am going to input: $\mu$, $\Sigma$, 1 (vector of ones), $H$, and $\alpha$, and try to get a value for $\gamma$ (also, the ' represents the transpose of the vector).

Equations:

$$ H = w(\gamma)'\mu + \Phi^{−1}(α)[w(\gamma)' \Sigma w(γ)]^{1/2}$$ $$ w(\gamma) = \frac{1}{\gamma}\Sigma^{-1}\Big[\mu − \Big(\frac{\mathbf{1}' \Sigma^{-1} \mu-\gamma}{\mathbf{1}' \Sigma^{-1} \mathbf{1}}\Big) \mathbf{1}\Big]$$

What I already tried: I am trying to solve this system using a python solver. But I am not able to do it since the matrix multiplication apparently becomes non-convex.

Question: How to solve this problem?

Obs: This question comes from the article Portfolio Optimization with mental accounts from Sanjiv Das, Harry Markowitz, Jonathan Scheid, and Meir Statman, Equations (7) and (8). The article is found at: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.410.8747&rep=rep1&type=pdf

Answer 1

I will use a case in which I randomly generate both $\Sigma$ and $\mu$ randomly in 4 dimensions

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.optimize import bisect

n = 4
np.random.seed(0)

S = np.random.uniform(0, 1, [n, n])
S = 0.5 * (S + S.transpose())
mu = np.random.uniform(-1, 1, n)

The trick with the transpose is just there to ensure that S is symmetric. I will choose $\alpha = 0.90$, and define the function $H$ as

Sinv = np.linalg.inv(S)
ones = np.ones_like(mu)

def H(gamma):

    # w
    w = (ones.dot(Sinv).dot(mu) - gamma) / (ones.dot(Sinv).dot(ones))
    w = Sinv.dot(mu - w) / gamma

    # H
    H = w.dot(mu) + norm.ppf(alpha) * np.sqrt(w.dot(S).dot(w))

    return H

This is a plot of $H = H(\gamma)$.

 g = np.linspace(4, 50, num = 500)
 h = np.array(map(H, g))
 i = h > 0
 plt.plot(g[i], h[i])
 plt.show()

The function is indeed sort of problematic for large values of $\gamma$ in the sense that any optimization algorithm (pretty much all of them) that uses derivatives will fail. But still you can use bisection to find the root

As an example I will use $H = 1.04$ (up in the flat region)

f = lambda x: H(np.exp(x)) - 1.04
r = bisect(f, 1.5, 6, xtol = 1e-12)
g0 = np.exp(r)
print g0, H(g0)

which generates the output

231.206991393 1.04

So the root is $\gamma = 231.20699$