portfolio optimization - markowitz

import numpy as np
from scipy import stats
import pandas as pd
from pandas_datareader import data as pdr
import datetime
import matplotlib.pyplot as plt
import fix_yahoo_finance as yf
yf.pdr_override()
import scipy.optimize as optimization


#Definições globais
janelaMovel = 1800
janelaMovelResultado = int (janelaMovel*0.69)

#esses são os ativos que escolhemos para nosso portifólio

stocks = ['PETR4.SA', 'ITUB4.SA', 'BBDC4.SA', 'BBAS3.SA', 'ABEV3.SA']


#nós usamos a base historica para aproximar média e variância: MPT depende da base histórica !!!
#fazer um parametro n para laço de repetição, mudar so n com intervalo que eu quero ~~ testar todos anos



#baixando a base de dados do Yahoo! Finance
def download_data(stocks, comeco, fim):
    data = pdr.get_data_yahoo(stocks, start=comeco, end=fim)['Adj Close']
    return data


#verificando se a distribuição de retornos é similar ao comportamento de uma normal
def show_data(data):
    data.plot(figsize=(10,5))
    daily_returns = (data/data.shift(1))-1
    daily_returns.hist(bins=100)
    plt.show()


#nos geralmente usamos logaritmo natural para normalizar os dados
def calculate_returns(data):
    returns = np.log(data/data.shift(1))
    return returns;

def plot_daily_returns(returns):
    returns.plot(figsize=(10,5))
    plt.show()

#printando a media e a covariancia dos ativos entre as datas;
def show_statistics(returns):
    print(returns.mean()*janelaMovelResultado)
    print(returns.cov()*janelaMovelResultado)

#definindo o peso dos ativos no portifolio de forma randomica
def initialize_weights():
    weights = np.random.random(len(stocks))
    weights /= np.sum(weights)
    return weights;

#retorno esperado com base nos pesos randomicos
def calculate_portfolio_return(returns, weights):
    portfolio_return = np.sum(returns.mean()*weights)*janelaMovelResultado
    print("Retorno esperado do portfolio:", portfolio_return)

#variancia esperada do mesmo portifolio randomico
def calculate_portfolio_variance(returns, weights):
    portfolio_variance = np.sqrt(abs(np.dot(weights.T, np.dot(returns.cov()*janelaMovelResultado,weights))))
    print("Variancia esperada do portfolio:", portfolio_variance)

def generate_portfolios(weights, returns):

    preturns = []
    pvariances = []

    #simulaçao de Monte-Carlo: no vamos gerar portifolios randomicos, com diferentes pesos de ativos
    for i in range(1000000):
        weights = np.random.random(len(stocks))
        weights/=np.sum(weights)
        preturns.append(np.sum(returns.mean()*weights)*janelaMovelResultado)
        pvariances.append(np.sqrt(np.dot(weights.T,np.dot(returns.cov()*janelaMovelResultado,weights))))

    preturns = np.array(preturns)
    pvariances = np.array(pvariances)
    return preturns,pvariances

def plot_portfolios(returns, variances):
    plt.figure(figsize=(10,6))
    plt.scatter(variances,returns,c=returns/variances,marker='o')
    plt.grid(True)
    plt.xlabel('Volatilidade Esperada')
    plt.ylabel('Retorno Esperado')
    plt.colorbar(label='Indice de Sharpe')
    plt.show()

# Ook, temos um resultado para simulaçao ... nos temos que encontrar agora o portifolio otimo
# Algumas tecnicas de otimização !!! scipy pode otimizar funções (minimo/maximo)
def statistics(weights, returns):
    portfolio_return=np.sum(returns.mean()*weights)*janelaMovelResultado
    portfolio_volatility=np.sqrt(np.dot(weights.T,np.dot(returns.cov()*janelaMovelResultado,weights)))
    return np.array([portfolio_return,portfolio_volatility,portfolio_return/portfolio_volatility])

# [2] media é o que queremos maximizar de acordo com o indice de sharpe
# note: maximizando f(x) é o mesmo que minimizar -f(x) !!!
def min_func_sharpe(weights,returns):
    return  -statistics(weights,returns)[2]

# what are the constraints? The sum of weights = 1 !!!  f(x)=0 this is the function to minimize
def optimize_portfolio(weights,returns):
    constraints = ( {'type':'eq','fun': lambda x: np.sum(x)-1})#, #a soma dos pesos é 1
                #{'type':'ineq','fun': lambda x: (10  - np.sum((x>=0)))},
            #{'type':'ineq','fun': lambda x: 8 - np.sum((x>=0))})

 #quantidade de ativos eh maior igual a 8
    bounds = tuple((0,1) for x in range(len(stocks))) #quando o peso for 1, estaremos investindo 100% em um único ativo
    optimum=optimization.minimize(fun=min_func_sharpe,x0=weights,args=returns,method='SLSQP',bounds=bounds,constraints=constraints)
    return optimum

# portifolio otimo de acordo com os pesos: 0 significa que n temos parcela alocada nesse ativo
def print_optimal_portfolio(optimum, returns):
    print("Pesos Otimos:", optimum['x'].round(5))
    print("Retorno esperado, volatilidade e indice de Sharpe:", statistics(optimum['x'].round(5),returns))

def print_result(optimum, returns):
    print("Pesos Otimos:", optimum['x'].round(5))
    print("Retorno, volatilidade e indice de Sharpe:", statistics(optimum['x'].round(5),returns))

def show_optimal_portfolio(optimum, returns, preturns, pvariances):
    plt.figure(figsize=(10,6))
    plt.scatter(pvariances,preturns,c=preturns/pvariances,marker='o')
    plt.grid(True)
    plt.xlabel('Volatilidade Esperada')
    plt.ylabel('Retorno Esperado')
    plt.colorbar(label='Indice de Sharpe')
    plt.plot(statistics(optimum['x'],returns)[1],statistics(optimum['x'],returns)[0],'g*',markersize=20.0)
    plt.show()

def future (optimum, stocks, inicio, fim):
    data = download_data(stocks, inicio, fim)
    returns2 = calculate_returns(data)
    #print_result(optimum, returns2)
    return statistics(optimum['x'].round(5), returns2)



if __name__ == "__main__":

    start = pd.to_datetime('2013-01-01')
    end = pd.to_datetime('2017-12-31')

    valoresAbsolutos = [1,0]
    n = 0



    while (start < (end - datetime.timedelta(days=janelaMovel))):
        base = pd.to_datetime(start - datetime.timedelta(days=janelaMovel))
        check = pd.to_datetime(start + datetime.timedelta(days=janelaMovel))

        print (base)
        print (check)

        print ("Realizando o download da base de dados")
        data = download_data(stocks, base, start)
        show_data(data)
        returns = calculate_returns(data)
        plot_daily_returns(returns)

        show_statistics(returns)
        weights=initialize_weights()
        calculate_portfolio_return(returns,weights)
        calculate_portfolio_variance(returns,weights)
        preturns,pvariances=generate_portfolios(weights, returns)
        plot_portfolios(preturns,pvariances)

        #print ('\n')
        print ("O portfolio otimo e: ")

        optimum=optimize_portfolio(weights,returns)
        print_optimal_portfolio(optimum, returns)
        show_optimal_portfolio(optimum, returns, preturns, pvariances)

        #print ('\n')
        #print ('\n')
        print ("O desempenho do portfolio otimo foi de: ")
        gain = future(optimum, stocks, start, check)



        n = n+1
        valoresAbsolutos[0] = valoresAbsolutos[0] + valoresAbsolutos[0]*gain[0]
        valoresAbsolutos[1] = float ((valoresAbsolutos[1] + gain[1]))


        start = check
        print ("----------------------------------------")

valoresAbsolutos[0] = valoresAbsolutos[0] - 1
valoresAbsolutos[1] = float ((valoresAbsolutos[1] + gain[1]))

print (valoresAbsolutos)