Fronteira Eficiente - Carteira de Markovitz

Resolveremos numericamente o problema para traçar a fronteira média-variância através da solução do problema de minimização da variância esperada da carteira, para um nível mínimo de retorno esperado como restrição.

Importação dos pacotes utilizados

Importamos os pacotes necessários e fazemos a inicialização da ferramenta que apresentará os gráficos.

In [1]:
import pandas_datareader.data as web
import matplotlib.pyplot as plt
import pandas as pd
import datetime as dt
import numpy as np
from scipy.optimize import minimize

from bokeh.io import show, output_notebook
from bokeh.plotting import figure
import bokeh.palettes as palettes
from bokeh.models import NumeralTickFormatter

output_notebook()
Loading BokehJS ...

Captura dos preços do ativos

Selecionamos 5 ativos negociados na bolsa brasileira para o problema. Capturamos os preços a partir do início do ano de 2016.

In [2]:
first_price_date = dt.datetime(2016, 1, 1)
last_price_date = dt.datetime(2016, 12, 7)
symbols = ['VALE5', 'PETR4', 'BVMF3', 'ITUB4', 'BBAS3']

data = pd.DataFrame()
for sym in symbols:
    data[sym] = web.DataReader(sym, 'google', first_price_date, last_price_date)['Close']
data = data.dropna()

days = data.index
print("From %s to %s \n" % (days[0].strftime("%d/%b/%y"), days[-1].strftime("%d/%b/%y")))
From 04/Jan/16 to 07/Dec/16 

Cálculo dos retornos diários

Note que o importante não é o preço dos ativos, e sim seus retornos.

In [3]:
daily_returns = np.log(data / data.shift(1))

Cálculo da matriz de covariância e dos retornos anualizados

Devemos sempre anualizar os retornos e a matriz de covariância. Para isso utilizamos a convenção de 252 dias úteis.

In [4]:
_252_DAYS = 252
covariance_matrix = np.asarray(_252_DAYS * daily_returns.cov())
annualized_returns = np.asarray(_252_DAYS * daily_returns.mean())

Matriz de covariância anualizada

In [5]:
print(pd.DataFrame(covariance_matrix, index=symbols, columns=symbols))
          VALE5     PETR4     BVMF3     ITUB4     BBAS3
VALE5  0.423061  0.237970  0.086874  0.100836  0.138764
PETR4  0.237970  0.425339  0.144189  0.162488  0.266100
BVMF3  0.086874  0.144189  0.132996  0.083185  0.128917
ITUB4  0.100836  0.162488  0.083185  0.132272  0.160479
BBAS3  0.138764  0.266100  0.128917  0.160479  0.337085

Retorno anualizado

In [6]:
print(pd.DataFrame(annualized_returns, index=symbols, columns=['Return']))
         Return
VALE5  1.112155
PETR4  0.908760
BVMF3  0.478882
ITUB4  0.431669
BBAS3  0.694936

Definição dos métodos de cálculo de variância (risco) e retorno da carteira

Definimos as equações de variância e retorno da carteira. Lembramos que o risco é medido pelo desvio padrão, que é a raiz quadrada da variância.

In [7]:
def portfolio_variance(w):
    _w = np.asarray(w)
    return _w.dot(covariance_matrix).dot(_w)

def portfolio_return(w):
    _w = np.asarray(w)
    return _w.dot(annualized_returns)

Definição do gradiente da variância para facilitar a convergência na minimização

In [8]:
def portfolio_variance_gradient(w):
    _w = np.asarray(w)
    return 2 * _w.dot(covariance_matrix)

Definição do método de minimização da variância com restrição de retorno mínimo para a carteira

Aqui definimos o método que, dado um retorno mínimo, minimiza a variância da carteira e retorna o vetor com a proporção de cada ativo. Devemos chamar atenção para as restrições:

$\bullet$ Se não há venda a descoberto, w deve ser não negativo.

$\bullet$ A soma dos elementos de w é sempre igual a 1.

$\bullet$ O retorno da carteira é sempre maior ou igual ao retorno mínimo passado como parâmetro.

In [9]:
def minimize_portfolio_variance(min_return, short_sale_allowed=True):
    assets_number = len(symbols)

    bnds = None if short_sale_allowed \
        else [(0, None) for i in range(assets_number)] # If short sale not allowed, lower bound = 0

    initial_guess = [1 / assets_number for i in range(assets_number)]
    cons = ({'type': 'eq', 'fun': lambda w: w.sum() - 1},
            {'type': 'ineq', 'fun': lambda w: portfolio_return(w) - min_return})
    return minimize(portfolio_variance, initial_guess,
                    constraints=cons,
                    bounds=bnds,
                    options={'disp': False},
                    method='SLSQP',
                    jac=portfolio_variance_gradient)

Definição do método que calcula a fronteira eficiente

O método que calcula a fronteira eficiente executa o método definido anteriormente para diversos valores de retorno mínimo. Assim, para cada iteração, temos um ponto da fronteira desejada.

In [10]:
def calculate_mv_frontier(short_sale_allowed=True):
    frontier = pd.Series()
    weights = pd.DataFrame()
    for min_return in np.linspace(0, 1.2, num=2000):
        result = minimize_portfolio_variance(min_return, short_sale_allowed=short_sale_allowed)
        if not result.success:
            continue
        w_frontier = result.x
        port_sigma = np.sqrt(portfolio_variance(w_frontier))
        port_return = portfolio_return(w_frontier)
        frontier.set_value(port_sigma, port_return)

        for i, symbol in enumerate(symbols):
            weights.set_value(port_sigma, symbol, w_frontier[i])

    return frontier, weights

Cálculo das fronteiras e das respectivas proporções dos ativos

Agora podemos calcular as fronteiras desejadas: Uma será calculada sem vendas a descoberto e a outra com vendas a descoberto. Medimos o tempo de execução para avaliar a eficiência do algoritmo de otimização.

In [11]:
%time frontier_ss_false, weights_ss_false = calculate_mv_frontier(short_sale_allowed=False)
%time frontier_ss_true, weights_ss_true = calculate_mv_frontier(short_sale_allowed=True)
CPU times: user 7.18 s, sys: 84.1 ms, total: 7.27 s
Wall time: 7.7 s
CPU times: user 8.57 s, sys: 87.8 ms, total: 8.66 s
Wall time: 9.32 s

Apresentação dos gráficos

Fronteiras com e sem vendas a descoberto

É interessante notar que a fronteira que permite vendas a descoberto é maior, por possibilitar mais combinações de ativos. Porém, como esperado, na região em que a carteira ótima possui apenas proporções positivas, as duas fronteiras se confundem.

Outra característica importante que deve ser notada é que o retorno máximo da carteira quando não é permitido venda a descoberto é o maior retorno dentre os ativos, alcançado quando a carteira está totalmente alocada neste ativo (proporção de 100%). Quando a venda a descoberto é permitida, alcançamos esse mesmo retorno com um risco inferior e com outros ativos na carteira.

Dica: Ao clicar no gráfico, selecione a opção de Zoom para verificar os níveis de risco.

In [12]:
p1 = figure(width=900, height=600, title="Efficient Frontier", toolbar_location="above")
p1.line(frontier_ss_true.index, frontier_ss_true.values, legend="Short Sale Allowed", color="green", line_width=2)
p1.line(frontier_ss_false.index, frontier_ss_false.values, legend="No Short Sale Allowed", color="blue", line_width=2)

p1.legend.location = "top_left"
p1.xaxis.axis_label = 'Risk'
p1.xaxis.formatter = NumeralTickFormatter(format="0.%")
p1.yaxis.axis_label = 'Return'
p1.yaxis.formatter = NumeralTickFormatter(format="0.%")

show(p1)

Proporção de cada ativo, com vendas a descoberto

Por não haver restrição nas proporções, notamos nas curvas abaixo uma aparente continuidade.

In [13]:
mypalette = palettes.Spectral5[0:len(symbols)]

p2 = figure(width=900, height=600, title="Efficient Frontier Weights (Short Sale Allowed)", toolbar_location="above")
for i, symbol in enumerate(symbols):
    p2.line(weights_ss_true[symbol].index, weights_ss_true[symbol].values, legend=symbol, color=mypalette[i], line_width=2)

p2.legend.location = "bottom_left"
p2.xaxis.axis_label = 'Risk'
p2.xaxis.formatter = NumeralTickFormatter(format="0.%")
p2.yaxis.axis_label = 'Weight'
p2.yaxis.formatter = NumeralTickFormatter(format="0.%")

show(p2)

Proporção de cada ativo, sem vendas a descoberto

A continuidade observada no gráfico anterior não pode mais ser observada quando colocamos a restrição das vendas a descoberto.

In [14]:
p3 = figure(width=900, height=600, title="Efficient Frontier Weights (No Short Sale Allowed)", toolbar_location="above")
for i, symbol in enumerate(symbols):
    p3.line(weights_ss_false[symbol].index, weights_ss_false[symbol].values, legend=symbol, color=mypalette[i], line_width=2)

p3.legend.location = "top_left"
p3.xaxis.axis_label = 'Risk'
p3.xaxis.formatter = NumeralTickFormatter(format="0.%")
p3.yaxis.axis_label = 'Weight'
p3.yaxis.formatter = NumeralTickFormatter(format="0.%")

show(p3)