Consideraciones
Como es algo extensa la respuesta (sólo se admite 30mil caracteres), el código completo está al final de este post, por tanto, sólo durante el desarrollo de esta respuesta, comentaré las funciones y líneas principales con los resultados respectivos para que sea más fácil su lectura.
Datos iniciales
Dado que se tiene muchos activos (50) tomamos 5 al alzar para comentar que se quiere hacer con ello
df = pd.read_csv("stocks.csv", sep=",", engine="python")
#���Field 1 la columna tiene caracteres extranios
df.columns = ['Field1' if i == 0 else x for i, x in enumerate(df.columns)]
df.set_index(df.columns[0], inplace=True)
# Seleccionamos 5 activos al azar
activos = np.random.permutation(df.columns)[:5]
plt.figure(figsize=(8,6))
plt.title('Historico de Algunos Activos', color='gray')
for activo in activos:
plt.plot(df[activo].to_numpy(), label=activo)
plt.ylabel('Precio del Activo', color='gray')
plt.xlabel('Observaciones', color='gray')
plt.legend(loc='upper left')
Como se ve los activos están a diferentes escalas (precios), hay que re-escalarlos, aunque en el pregunta para ello se emplea returns = stock_prices.pct_change()
en esta situación en particular es aconsejable usar logaritmos... ¿Por qué usar rendimiento logarítmicos? con ello se obtiene la siguiente gráfica:
df = df / df.shift(1)
df.dropna(inplace=True)
log_df = np.log(df)
El Modelo de Frontera Eficiente de Markowitz, requiere una matriz de covarianza (C
) y unos rendimientos históricos (p
) que dependiendo de su periodo (diarios semanales, mensual, etc.) estos pueden ser anualizados... en este caso concreto 256 días por año
def metricas_historicas_portafolio(portafolio, dias_anual):
"... resto del codigo"
return p, C, numero_de_activos
p, C, numero_de_activos = metricas_historicas_portafolio(log_df, 252)
Rendimiento y Volatilidad
De acuerdo a las colocaciones de nuestro fondo de inversión en un determinado portafolio se obtendrá un rendimiento (mu) y volatilidad (riesgo o sigma) del mismo. Dichas colocaciones corresponderán a unos pesos de relevancia que le hemos dados a cada uno de los activos en nuestro portafolio, estos vienen dados por las siguientes expresiones de matrices:
Al usar numpy
dichas expresiones se pueden implementar tal cual, pero con una pequeña variación en el orden, por la forma en como fue implementada la función metricas_historicas_portafolio
... más detalle en este post en inglés:
def resultados_portafolio(p,w,C):
mu = w * p.T # Rendimiento Esperado
sigma = np.sqrt(w * C * w.T) # Volatilidad
return mu, sigma
Ratio de Sharpe
Es una medida del exceso de rendimiento por unidad de riesgo de una inversión wikipedia
Es útil porque podemos obtener el desempeño de nuestro portafolio, es decir a más grande mejor.
Debido a su importancia hay una función implementada para obtener métricas sobre un conjunto de portafolios
def rsharpe_maximo(sharper_ratios, rendimientos, volatilidades):
"... resto del codigo"
return maximo_sharpe_ratio, maximo_sharpe_ratio_volatilidad, maximo_sharpe_ratio_rendimiento
Montecarlo
Supongamos que el inversionista quiere colocar USD $100
este portafolio con 50 activos
... una posible e "ingenua" combinación sería colocar USD $2 en cada activo
... es decir los pesos sería [0.02, 0.02, 0.02, ..., 0.2] para cada activo
Generalmente se simulan los pesos con una distribución uniforme o distribución normal estándar para obtener pesos "aleatorios" y ver posibles portafolios...
Lo anterior no está mal si son pocos activos... de lo contrario se está obviando la ley débil de los números grandes es decir se estaría emulando el comportamiento de un inversionista ingenuo
Veamos la siguiente función implementada para explicar de que trata
def simular_pesos_portafolio(numero_de_activos):
pesos = np.random.random(numero_de_activos)
pesos *= sum([np.random.binomial(1, 0.08, numero_de_activos) for _ in range(2)])
pesos = np.asmatrix(pesos / sum(pesos) )
return pesos
Para evitar lo anterior se añade un generador binomial de 0s
y 1s
con una probabilidad de 0.08 a efecto de que no se invierta en algunos activos... que sucede si todos salen 0s
... como se verá más adelante el Sharpe Ratio
será nan
... hay que preocuparse... la respuesta corta y directa es no
La idea de simular para este caso concreto, es ver el comportamiento de las posibles opciones que tenemos en la mesa, ya que al existir un modelo matemático robusto para optimizar nuestro portafolio, este debe ser usado por su eficiencia en computo y exactitud.
Simulamos 1000
portafolios con la función simulacion_de_portafolios
que llama en cada iteración a la función simular_portafolio
que su vez depende de la función simular_pesos_portafolio
pesos, rendimientos, volatilidades, sharper_ratios = simulacion_de_portafolios(
numero_de_portafolios=1000,
p=p,
C=C,
numero_de_activos=numero_de_activos,
libre_riesgo=0
)
También se obtienen las métricas del ratio de sharpe
maximo_sharpe_ratio, maximo_sharpe_ratio_volatilidad, maximo_sharpe_ratio_rendimiento = \
rsharpe_maximo(sharper_ratios, rendimientos, volatilidades)
Si se cambia la semilla np.random.seed(9062020)
también cambiará las métricas... de lo que se deduce cuál es el "óptimo"
Solvers ("Solucionadores")
Si bien al leer u oír solver se le asocia a Excel, hay que tener presente que es una tecnología que permite minimizar o maximizar una función objetivo para un conjunto de restricciones dadas.
En el enlace proporcionado se da una lista de "Solvers" no lineales... sobre el particular voy a ser honesto, no les he probado todos, por lo general trabajo con OR Tools for Python (Herramientas de Investigación de Operaciones para Python) de Google, pero esta tecnología que pertenece a un conjunto de herramientas de Inteligencia Artificial de Google, que es Open Source, a la fecha sólo está orientada a modelos lineales.
Para esta respuesta estoy empleando el solver cvxopt, porque hay varios ejemplos proporcionados en la red sobre este caso concreto, incluso hay un código ejemplo para la frontera eficiente, el cuál he tomado pero haciéndole algunas modificaciones, para emplear la función resultados_portafolio
implementada previamente
def portafolio_optimo(p, C, numero_de_activos):
# Se establece saltos discretos para hallar la
# la frontera eficiente estos seran los
# `targets` u objetivos que se fijan para optimizar
N = 100
n = numero_de_activos
#mus = np.power(10, 5 * np.arange(N) / N - 1) # tiene que ser lista .tolist()
mus = [ 10**(5.0*t/N-1.0) for t in range(N) ]
# convertir el p y C a matrices del tipo cvxopt
# en el caso de p se trabaja con su transpuesta
pbar = opt.matrix(p.T)
S = opt.matrix(C)
# Crear las matrices de restricciones
# Gx <= h
G = -opt.matrix(np.eye(n)) # matriz identidad negativa n x n
h = opt.matrix(0.0, (n ,1))
# Ax = b
A = opt.matrix(1.0, (1, n))
b = opt.matrix(1.0) # La suma de los pesos es 1
# Calcular los pesos de la frontera eficiente
# Empleando Programacion Cuadratica
# Pero primero silenciamos el solver (es opcional)
solvers.options['show_progress'] = False
portafolios = [solvers.qp(mu*S, -pbar, G, h, A, b)['x'] for mu in mus]
# Calcular los rendimientos y volatilidades o riesgos
# para la frontera eficiente
# Estas implementaciones funcionan... hay que importar "blas"
# pero la implementacion que esta fuera de este codigo
# requiere que se redimensionen los rendimientos, volatilidades
# para que funcionen
#rendimientos = [blas.dot(pbar, x) for x in portafolios]
#volatilidades = [np.sqrt(blas.dot(x, S*x)) for x in portafolios]
rendimientos, volatilidades = zip(*[
resultados_portafolio(p, np.array(w).T, C) for w in portafolios
])
# Calcular el portafolio optimo
#pesos = solvers.qp(opt.matrix(x1 * S), -pbar, G, h, A, b)['x'] # No funciona para lo que se quiere hacer mas adelante
pesos = [np.asarray(x) for x in portafolios]
pesos = np.asarray(pesos)
rendimientos = np.asarray(rendimientos)
volatilidades = np.asarray(volatilidades)
return pesos, rendimientos, volatilidades
w_optimos, mu_optimos, sigma_optimos = portafolio_optimo(p, C, numero_de_activos)
La implementación anterior corresponde a un desarrollo formal de
a:
Con ellos se genera los puntos para la Frontera Eficiente con los respectivos pesos que dan dichas coordenadas (cada punto es un portafolio)
Pero antes de ellos quiero comentar la siguiente linea
mus = [ 10**(5.0*t/N-1.0) for t in range(N) ]
Esto contiene los rendimientos objetivos o targets los cuales crecen de forma exponencial... lo cuál es válido para un inversionista que busca obtener un desempeño alto de su dinero (similar a las startups de TI que crecen x2 x3 x4 anualmente... los llamados unicornios)
Y ¿por qué usar "targets"?... porque al ser la Frontera Eficiente una curva en ella hay infinitos puntos o infinitos portafolios que cae en dicha curva.
En tal sentido, no hay nada esotérico con esos números por defectos, puedes cambiarlos para obtener una curva más fina si así lo deseas (N>100) saltos discretos más pequeños, etc, pero eso sí... que crezca de forma exponencial.
Pesos (Colocación)
Cuando se simula hay rendimientos positivos y negativos para un portafolio cualquiera. En el caso de la Frontera Eficiente se aprecia que todos los rendimientos...
Pero esto implica ¿qué los pesos o colocaciones sean necesariamente positivas?
w_optimos, mu_optimos, sigma_optimos = portafolio_optimo(p, C, numero_de_activos)
print("\nVerificar cuales suman 1: ")
print(np.array([x.sum() for x in w_optimos]))
# [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
# 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
# 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
# 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
# 1. 1. 1. 1.]
filtrar_pesos_positivos = np.array([(x>=0).all() for x in w_optimos])
print("\nVerificar que todos los pesos sean >= 0: ")
print(filtrar_pesos_positivos)
# [ True True True True True True True True True True True True
# True True True True True True True True True True True True
# True True True True True True True True True True True True
# True True True True True True True True True True True True
# True True True True True True True True True True True True
# True True True True True True True True True True True True
# True True False False False False True False False False False True
# True True False False False False False True True True True False
# False False True True]
El hecho de tener pesos (colocaciones) negativas ya fue ampliamente explicado acá
Solve(covMat) retorna que system is computationally singular: reciprocal condition number
No obstante, se puede descartar dichos portafolios ya que fueron generados 100
quedándonos sólo 84
con esta sentencia:
filtrar_pesos_positivos = np.array([(x>=0).all() for x in w_optimos])
y actualizando los puntos de la Frontera Eficiente generados previamente
w_optimos = w_optimos[filtrar_pesos_positivos]
mu_optimos = mu_optimos [filtrar_pesos_positivos]
sigma_optimos = sigma_optimos[filtrar_pesos_positivos]
Ratio de Sharpe - de nuevo ;) -
Con la función rsharpe_maximo
y habiéndose descartado aquellos portafolios con pesos negativos, se busca le máximo Ratio de Sharpe para la Frontera Eficiente
maximo_rsharpe, maximo_rsharpe_volatilidad, maximo_rsharpe_rendimiento = \
rsharpe_maximo(rsharpe_optimos, mu_optimos, sigma_optimos)
print(maximo_rsharpe)
# 1.27
Resultado Final (Todo Junto)
Como se sobre-escribió las variables del portafolio_optimo
se ejecuta (sé que se puede pulir esta parte para reutilizar el código pero la idea en este caso es mostrar un paso a paso para luego optimizar la escritura del código si se requiere)
w_optimos, mu_optimos, sigma_optimos = portafolio_optimo(p, C, numero_de_activos)
libre_riesgo = 0
rsharpe_optimos = (mu_optimos - libre_riesgo) / sigma_optimos
rsharpe_optimos = rsharpe_optimos.reshape((1,-1)).reshape((1,-1)) # quitarle dimensiones
Como los resultados no van a cambiar (a diferencia de un simulación) y ya se tiene el ratio de sharpe máximo para aquellos portafolios cuyos pesos son positivos, se emplea el siguiente código para capturar su posición en el gráfico final. Aquí el fragmento del código relevante:
# Optimo Sharpe Ratio Frontera Eficiente
idx_rshape_optimo = np.where(rsharpe_optimos[0] == maximo_rsharpe)
plt.scatter(
sigma_optimos.reshape((1, -1))[0][idx_rshape_optimo], # eje volatilidad
mu_optimos.reshape((1,-1))[0][idx_rshape_optimo], # eje rendimientos
c='#90ff1e', s=100, edgecolors='gray', label=f'Sharpe Ratio Optimo = {maximo_rsharpe:.4f}'
)
Notas
Todos los gráficos (excepto los dos iniciales) está anualizados
(rendimiento y volatilidad para 252 días al año)
Se ha empleado reshape((1,-1))
en las variables w_optimos, mu_optimos, sigma_optimos
esto por que al trabajar con la función
resultados_portafolio
y con algebra lineal se obtienen en
dimensiones que matplotlib
no puede leer directamente. Sí en la
función portafolio_optimo
quisieras implementar las funciones tal
cual está en la documentación de cvxopt:
rendimientos = [blas.dot(pbar, x) for x in portafolios]
volatilidades = [np.sqrt(blas.dot(x, S*x)) for x in portafolios]
Vas a tener que tener en cuenta que no se obtendrán las mismas
dimensiones para mu_optimos, sigma_optimos
con lo cual tendrás que
re-implementar toda la parte de abajo (después de la función
portafolio_optimo
) para que calcen las dimensiones de los arreglos.
Debido a que el archivo contiene el siguiente carácter ���Field 1
no sé en que versión lo estás ejecutando quizás esto pueda ayudar ¿Qué es el error “SyntaxError: Non-ASCII character '\xc2', but no encoding declared”?
o en el terminal ejecutar la siguiente instrucción con el nombre de tu fichero
sin la extensión .py
que contiene el código más abajo posteado en mi caso el archivo se llama modeloMarkowitz.py:
$ PYTHONIOENCODING=utf8 python3 -c "import modeloMarkowitz"
Código
Como el post es extenso (máximo 30mil carácteres), el código completo que desarrollé, está en mi repositorio de GitHub ModeloMarkowitz