# -*- coding: utf-8 -*-
"""
Portfolio : random and optimal
"""
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import softmax
from tqdm import tqdm
# version avec csv sur le site du cours
data = pd.read_csv('close_cac40_historical.csv',sep=';', index_col = 0)
data.head()
data.keys()
data0=data.copy()

#select first nb columns
nb=40
data = data.iloc[:,0:nb]

returns= np.log(data/data.shift(1))#rdt log
#calcul de la moyenne des returns, et la covariance
mean_returns=returns.mean()*250#annualized
cov_matrix= returns.cov()*250#annualized
nb = mean_returns.shape[0]

nbportef=3000
rdt,vol=[],[]
for ik in tqdm(range(nbportef)):#use tqdm for progression bar if nbportef = high
#for ik in range(nbportef):
    #select portfolios at random and plot them
    allocation=softmax(np.random.randn(nb))# pour rendre la somme=1 : softmax
    rdt.append(allocation@mean_returns)
    vol.append(np.sqrt(allocation@cov_matrix@allocation))    

plt.plot(vol,rdt,"*")

#compute a,b, inverse_cov
inverse_cov= np.linalg.inv(cov_matrix)
a=np.ones(nb)@inverse_cov@np.ones(nb)
b= mean_returns@inverse_cov@np.ones(nb)
tmpvect= mean_returns-(b/a)*np.ones(nb)
normvect= np.sqrt(tmpvect@inverse_cov@tmpvect)

volrange=np.linspace( 1/np.sqrt(a)+1.e-6,1.1*np.max(vol))
optimal_returns= b/a + np.sqrt(volrange**2-1/a)*normvect
plt.plot(volrange,optimal_returns,"-g")
plt.xlabel("volatility")
plt.ylabel("mean return")