# -*- coding: utf-8 -*-
"""
1/ Telechargement cours composantes CAC40

2/ historique des prix et des rendements

3/ test de normalité avec scipy.stats.normaltest

"""


import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

data_p=pd.read_pickle("close_cac40_historical.pkl")  
data_p.head()
data_p.keys()

# version avec csv sur le site du cours
data = pd.read_csv('close_cac40_historical.csv',sep=';', index_col = 0)
#calcul des rendements
rendements = np.log(data/data.shift(1))
#annuel
rendements_an = np.log(data/data.shift(255))

dim=39
rendements_dim=rendements_an.iloc[:,0:dim]

M=rendements_dim.mean()
Sigma=rendements_dim.cov()
SigmaInv = np.linalg.inv(Sigma)

a = np.sum(SigmaInv)
b = np.ones(dim)@SigmaInv@M
print("a=",a,' b=',b)

xrange=np.linspace(1/np.sqrt(a)+1.e-6,
                   4/np.sqrt(a),50)

factor= M-b/a
norm_factor= np.sqrt(factor@SigmaInv@factor)
optimal=b/a+norm_factor*np.sqrt(xrange**2-1/a)

plt.figure('portef')
plt.plot(xrange,optimal)


#calcul des portef au hasard
no_portef=1000
return_portef=[]
volatility_portef=[]
for _ in range(no_portef):
    #choix loi uniforme
    weights= np.random.uniform(-4./dim,4./dim,dim)
    weights = weights - np.mean(weights) + 1.0/dim

    #choix loi normale

    #calculer return, volatility
    return_portef.append(weights @ M)
    
    current_volatility= np.sqrt(weights.T @ Sigma @ weights)
    volatility_portef.append(current_volatility)
    
#plot
plt.figure('portef')
plt.plot(volatility_portef,return_portef,'.b',
         xrange,optimal,'-*g')
plt.xlabel('sigma')
plt.ylabel('return')
plt.title('Return vs. sigma for n='+np.str(dim))
plt.show()

#alternatives np.random.dirichlet mais sans decouvert