# -*- coding: utf-8 -*-
"""
A FAIRE:
    
    1/ telecharger les cours des composantes du CAC40 (ou DJIA, ou ...)
    en prix journaliers (au moins le "close") si possible depuis 5 ans
    
    Idee: utiliser yahoo comme source et par exemple yfinance 
    "pip install yfinance"

    Sinon prendre le fichier 
    https://turinici.com/wp-content/uploads/cours/common/close_cac40_historical.csv
    
    1'/ ordonner par date croissante
    
    2/ dessiner l'histograme du prix et du rendement (avec "log" et "actuariel")
    
    3/ tester la normalité des prix, celle des rendements log 
    et "actuariel" (par exemple utiliser scipy.stats.normaltest)
"""


import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from  scipy.stats import normaltest,kstest
from scipy.special import softmax


# take csv from www cours
data = pd.read_csv('close_cac40_historical.csv',sep=';', 
                   index_col = 'Date')
data.head()
data.sort_values(by='Date',inplace=True)#alternative "sort_index"
data.head()
data.tail()

# plot histogram of prices
_=data.hist(bins=30, figsize = (15,15))

def normality_test(data,kolmogorov_smirnov=False,level=0.01,print_results=True):
    """
    Tests normality of each column of dataframe "data". 
    Inputs:    
        kolmogorov_smirnov= false: use "normaltest", otherswise use kstest, both from scipy.stats
        level = p-value threshold level for the conclusions
        
    Outputs:  the number of yes/no in the results
    """

    pvalues=[]
    for cols in data.keys():
        # normality test 
        if(kolmogorov_smirnov):
            pv=kstest(data[cols].dropna(),'norm').pvalue
        else:
            pv=normaltest(data[cols],nan_policy='omit').pvalue
            
        pvalues.append(pv)
        if(pv<level):
            res="not normal"
        else:
            res='normal'
        print("Test pval=",pv, 'res=',res)
    normalok=sum([1 for pv in pvalues if pv >=level])    
    normalnotok=sum([1 for pv in pvalues if pv <level])    
    if(print_results):
        print("no. of normal = ",normalok)    
        print("no. of not normal = ",normalnotok)
    return normalok,normalnotok

normality_test(data)

#compute returns, plot hist, test
#returns= (data.diff()/data.shift(1)).shift(-1)# actuarial return
returns= np.log(data/data.shift(1))#rdt log
_=returns.hist(bins=int(np.sqrt(returns.shape[0])), figsize = (15,15))

normality_test(returns.tail(25*3))#test last 3 months

###########################################################
print('normality tests for increments, not returns!!')
increments= data-data.shift(1)
_=increments.hist(bins=int(np.sqrt(increments.shape[0])), figsize = (15,15))

normality_test(increments.tail(25*3))
########################################################################



#%%
nb=10  #will work with 10 stocks
all_returns=returns.copy()#backup
nb_all=all_returns.shape[1]
assert nb<=nb_all, "wrong number of stocks"
#returns_small=all_returns.iloc[:,np.random.choice(nb_all,nb)]
returns_small=all_returns.loc[:,
   np.random.choice(all_returns.keys(),nb,replace=False)]

#%%
#####################################
#calcul de la moyenne des returns, et la covariance
mean_returns=returns_small.mean()
cov_matrix= returns_small.cov()
plt.figure('portef')
rdt_list=[]
std_list=[]
for _ in range(500):
#    allocation=np.random.randn(nb)#a rendre la somme=1 ex: softmax
#    allocation=allocation/np.sum(allocation)
    allocation=softmax(np.random.randn(nb))# pour rendre la somme=1 : softmax
    rdt_port = allocation@mean_returns
    std_port = np.sqrt(allocation@cov_matrix@allocation)
    rdt_list.append(rdt_port)
    std_list.append(std_port)


inverse_cov = np.linalg.inv(cov_matrix)
# compute and draw the efficient frontier on the same graph
onesM=np.ones_like(mean_returns)
a=np.float(onesM @ inverse_cov @ onesM)
b=np.float(onesM @ inverse_cov @ mean_returns)

# dessiner la frontiere et sa symetrique par rapport a l'origine
sigmarange= np.linspace(1./np.sqrt(a)+1.e-10, 1.1*np.max(std_list),47)
#compute the optimal portfolio return for these values 
factor = np.sqrt((mean_returns- b/a*onesM) @ inverse_cov @ (mean_returns- b/a*onesM))
optimal_return = b/a + np.sqrt(sigmarange**2-1./a)*factor



plt.scatter(std_list,rdt_list)
plt.plot(sigmarange,optimal_return,'r-')
plt.xlabel('std')
plt.ylabel('rdt')
#plt.xlim([0,.2])
#plt.ylim([-.05,.05])