Resolución de un sistema de dos ecuaciones con multiplicación matricial no convexa para la optimización de MV

Question

Escenario: Estoy tratando de hacer una variación de la optimización de MV para una cartera. En este caso, ya tengo un vector de rendimientos medios ( $\mu$ ), un vector de unos, una matriz de covarianza ( $\Sigma$ ), y $\phi^{-1}$ que es la inversa de la función de distribución acumulativa normal (función de probabilidad alfa).

El problema: A partir de las dos ecuaciones indicadas a continuación, voy a introducir: $\mu$ , $\Sigma$ 1 (vector de unos), $H$ y $\alpha$ y tratar de obtener un valor para $\gamma$ (además, el ' representa la transposición del vector).

Ecuaciones:

$$ 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]$$

Lo que ya he probado: Estoy tratando de resolver este sistema utilizando un solucionador de python. Pero no soy capaz de hacerlo ya que la multiplicación de la matriz aparentemente se convierte en no convexa.

Pregunta: ¿Cómo resolver este problema?

Obs: Esta pregunta proviene del artículo Portfolio Optimization with mental accounts from Sanjiv Das, Harry Markowitz, Jonathan Scheid, and Meir Statman, Equations (7) and (8). El artículo se encuentra en: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.410.8747&rep=rep1&type=pdf

Answer 1

Utilizaré un caso en el que genero aleatoriamente ambos $\Sigma$ y $\mu$ al azar en 4 dimensiones

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)

El truco de la transposición es sólo para asegurar que S es simétrica. Elegiré $\alpha = 0.90$ y definir la función $H$ como

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

Esta es una parcela de $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()

La función es, en efecto, algo problemática para valores grandes de $\gamma$ en el sentido de que cualquier algoritmo de optimización (prácticamente todos) que utilice derivadas fallará. Pero aún así se puede utilizar bisección para encontrar root

Como ejemplo utilizaré $H = 1.04$ (arriba en la región plana)

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)

que genera la salida

231.206991393 1.04

Así que root es $\gamma = 231.20699$