Home / Questions / 361549
Accepted
Revolucion for Monica
Asked: 2020-06-03 01:00:16 +0800 CST 2020-06-03 01:00:16 +0800 CST

Maximize a variable with the constraint of a fixed volatility

  • 772

I am trying to create an Efficient Frontier based simulated portfolio optimization on 50 stocks, which you can find in csv here . Now I know how to get the best portfolio but it keeps the free weights and I'm looking for an optimization that keeps this variable positive.

I get my efficient frontier by getting weights with a fixed network of volatiles σ p 1 ,... σ p n . So for each σ p i , I maximize the expected returns with the constraint that the volatility is not greater than σ p i , to obtain μ p i . Then σ p i , μ p i are n points on the efficient frontier.

enter image description here

In other words, we can maximize the returns r that are equal to the weights of the individual results of each action RW . This leads to the following optimization problem (I reduced it to two variables for simplicity):

enter image description here

Which, in matrix form, should be more or less:

enter image description here

Where COV is the covariance matrix between all assets.

However, I was only able to program it with no weight restrictions W .

Here, the minimization of a quadratic objective under linear constraints can be obtained by solving a linear system. Which is really easy to calculate and not that hard to understand. To get the best acuity ratio give the following:

def efficient_portfolios(returns, risk_free_rate, sigma, mu, e):
    weights_record = []
    volatilities = []
    results = np.zeros((3,len(returns)))
    i = 0
    for portfolio_return in returns:
        A = np.block([[2*sigma, mu, e], [mu.T, 0, 0], [e.T, 0, 0]])
        b = np.zeros(n+2)
        b[n] = portfolio_return
        b[n+1] = 1
        w = np.linalg.solve(A, b)[:n]
        weights_record.append(w)
        portfolio_std_dev = np.sqrt( w.T @ sigma @ w )
        volatilities.append(portfolio_std_dev)
        results[0,i] = portfolio_std_dev
        results[1,i] = portfolio_return
        results[2,i] = (portfolio_return - risk_free_rate) / portfolio_std_dev
        i+=1
    return results, weights_record, volatilities

def display_simulated_ef_with_random(mean_returns, risk_free_rate, sigma, mu, e, df):
    results, weights, volatilities = efficient_portfolios(mean_returns,risk_free_rate, sigma, mu, e)
    max_sharpe_idx = np.argmax(results[2])
    sdp, rp = results[0,max_sharpe_idx], results[1,max_sharpe_idx]
    max_sharpe_allocation = pd.DataFrame(weights[max_sharpe_idx],index=df.columns,columns=['allocation'])
    max_sharpe_allocation.allocation = [round(i*100,2)for i in max_sharpe_allocation.allocation]
    max_sharpe_allocation = max_sharpe_allocation.T

    min_vol_idx = np.argmin(results[0])
    sdp_min, rp_min = results[0,min_vol_idx], results[1,min_vol_idx]
    min_vol_allocation = pd.DataFrame(weights[min_vol_idx],index=df.columns,columns=['allocation'])
    min_vol_allocation.allocation = [round(i*100,2)for i in min_vol_allocation.allocation]
    min_vol_allocation = min_vol_allocation.T

    print("-"*80)
    print("Maximum Sharpe Ratio Portfolio Allocation\n")
    print("Annualised Return:", round(rp,2))
    print("Annualised Volatility:", round(sdp,2))
    print("\n")
    print(max_sharpe_allocation)
    print("-"*80)
    print("Minimum Volatility Portfolio Allocation\n")
    print("Annualised Return:", round(rp_min,2))
    print("Annualised Volatility:", round(sdp_min,2))
    print("\n")
    print(min_vol_allocation)

    plt.figure(figsize=(10, 7))
    plt.scatter(results[0,:],results[1,:],c=results[2,:],cmap='YlGnBu', marker='o', s=10, alpha=0.3)
    plt.colorbar()
    plt.scatter(sdp,rp,marker='*',color='r',s=500, label='Maximum Sharpe ratio')
    plt.scatter(sdp_min,rp_min,marker='*',color='g',s=500, label='Minimum volatility')
    plt.title('Simulated Portfolio Optimization based on Efficient Frontier')
    plt.xlabel('annualised volatility')
    plt.ylabel('annualised returns')
    plt.legend(labelspacing=0.8)

    return max_sharpe_allocation, min_vol_allocation

stock_prices = pd.read_csv('stocks.csv', index_col=0)
returns = stock_prices.pct_change()
mu = 252 * returns.mean().values
sigma = 252 * returns.cov().values
n = mu.shape[0]

# add risk free asset to mu/sigma
risk_free_rate = 0.0178
z = np.zeros((n,1))
#mu = np.block([mu, risk_free_rate])
#sigma = np.block([[sigma, z], [z.T, 0]])
#n = mu.shape[0]

# solve minimize w'∑w subject to μ'w = r, e'w=1 for varying r
mu = np.expand_dims(mu, axis=1)
e = np.ones((n,1))
returns = np.linspace(risk_free_rate, np.max(mu))


max_sharpe_al, min_vol_al = display_simulated_ef_with_random(returns, risk_free_rate, sigma, mu, e, stock_prices)

And gives the following plot and portfolio:

--------------------------------------------------------------------------------
Maximum Sharpe Ratio Portfolio Allocation

Annualised Return: 0.63
Annualised Volatility: 0.23


               DD  ADBE  ATVI   APD    NVS      A   ADI    AVB    AYI   AAN  \
allocation -19.33  0.03 -0.32  29.3  12.65 -14.57  2.85 -25.28 -13.17  2.77   

            ...   SWKS    NOV  KMT   MDT   RIO   PSA   STE  POWI  VALE   TX  
allocation  ... -15.61 -10.08 -7.2 -3.16  7.57 -9.39  7.93  5.13  1.07  8.4  

[1 rows x 51 columns]
--------------------------------------------------------------------------------
Minimum Volatility Portfolio Allocation

Annualised Return: 0.03
Annualised Volatility: 0.13


             DD  ADBE  ATVI   APD   NVS     A   ADI   AVB   AYI   AAN  ...  \
allocation -0.6 -7.11  5.36  3.81  22.9 -3.69  7.37 -1.27 -1.13 -0.16  ...   

            SWKS   NOV   KMT   MDT   RIO    PSA   STE  POWI  VALE    TX  
allocation  -6.4 -0.25 -9.24  6.15  4.41  19.86 -1.31 -0.23 -2.99  6.05

enter image description here

But how would you do if we needed to always be positive? I know that the solution cannot be obtained simply by solving a linear system. You would need to use a quadratic optimization solver .

python

1 Answers

  1. Best Answer

    Considerations

    As the answer is somewhat long (only 30 thousand characters are allowed), the complete code is at the end of this post, therefore, only during the development of this answer, I will comment on the main functions and lines with the respective results to make it easier your reading.

    Initial data

    Since there are many assets (50) we take 5 randomly to comment on what we want to do with it

    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')

    enter image description here

    As you can see, the assets are at different scales (prices), they must be re-scaled, although in the question it is used returns = stock_prices.pct_change()in this particular situation, it is advisable to use logarithms... Why use logarithmic returns? With this, the following graph is obtained:

    df = df / df.shift(1) 
df.dropna(inplace=True)
log_df = np.log(df)

    enter image description here

    The Markowitz Efficient Frontier Model requires a covariance matrix ( C) and historical returns ( p) that depending on their period (daily, weekly, monthly, etc.) these can be annualized... in this specific case256 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)

    Yield and Volatility

    According to the placements of our investment fund in a certain portfolio, a return (mu) and volatility (risk or sigma) of the same will be obtained. These placements will correspond to relevant weights that we have given to each of the assets in our portfolio, these are given by the following matrix expressions:

    enter image description here enter image description here

    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

    enter image description here

    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)

    enter image description here

    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

    enter image description here enter image description here

    a:

    enter image description here

    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.

    enter image description here

    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

    enter image description here

    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}'
)

    enter image description here

    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”?

    • or in the terminal execute the following instruction with nombre de tu ficherowithout the extension .pythat contains the code below posted in my case the file is called ModeloMarkowitz.py:

      $ PYTHONIOENCODING=utf8 python3 -c "import modeloMarkowitz"

    Code

    As the post is long (maximum 30 thousand characters), the complete code that I developed is in my GitHub repository ModeloMarkowitz

    • 6