#!/usr/bin/env python
# coding: utf-8

# # 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()


# ## 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")))


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


# ### Retorno anualizado

# In[6]:


print(pd.DataFrame(annualized_returns, index=symbols, columns=['Return']))


# ## 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]:


get_ipython().run_line_magic('time', 'frontier_ss_false, weights_ss_false = calculate_mv_frontier(short_sale_allowed=False)')
get_ipython().run_line_magic('time', 'frontier_ss_true, weights_ss_true = calculate_mv_frontier(short_sale_allowed=True)')


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