"""
A FAIRE:
    
    1/ download quotations of CAC40 components (or DJIA, or ...) in daily
    prices (at least the "close") if possible at lest 5 years
    
    Idea: use yahoo e.g., yfinance package 
    "pip install yfinance"

    Otherwise take the file 
    https://turinici.com/wp-content/uploads/cours/common/close_cac40_historical.csv
    
    1'/ order by increasing date
    
    2/ plot price histogram and returns  (with "log" and/or "actuarial")
    
    3/ test normality of : prices, log returns, actuarial returns
    for instance can use scipy.stats.normaltest
    
    ATTENTION: replace "None" by your code !!
"""


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()
# order by increasing date, keep variable 'data'
data.sort_index(ascending=True,inplace=True)
#data.sort(by='Date',inplace=True)
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):
		# put in pv the p-value of the Kolmogorov-Smirnov test
            pv=kstest(data[cols].dropna(),'norm').pvalue
        else:
		#  put in pv the p-value of the default scipy.stats normality test
            pv=normaltest(data[cols].dropna()).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)

# use 'data' to compute returns
# returns= data.pct_change() #actuarial
returns= np.log(data/data.shift(1))

_=returns.hist(bins=int(np.sqrt(returns.shape[0])), figsize = (15,15))

normality_test(returns.tail(25*3))#test last 3 months
#normality_test(returns)#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.loc[:,
   np.random.choice(all_returns.keys(),nb,replace=False)]

#%%
#####################################
#compute avg and cov of returns
mean_returns=returns_small.mean()
cov_matrix= returns_small.cov()
plt.figure('portef')
rdt_list=[]
std_list=[]
for _ in range(500):
#A FAIRE : sample at random some "allocation". Use your own guess, only constraints is
#to have sum=1
# suggestions:  normal r.v. with a normalization factor for the norm;  softmax, ... etc
# "allocation"  format : numpy.ndarray of shape=(nb,)
    allocation=np.random.randn(nb)
    allocation=allocation/np.sum(allocation)
#    allocation=softmax(np.random.randn(nb))
#     allocation = np.random.dirichlet(onesM,size=1)[0]
    
    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 FAIRE : calculer 'a' et 'b' comme dans les formules de cours
a=None
b=None

# plot the frontier and its symmetric w/r to origin
sigmarange= np.linspace(1./np.sqrt(a)+1.e-10, 1.1*np.max(std_list),47)
#TODO : compute the return of the optimal portfolio for sigma in sigmarange
# will use the "factor" auxiliary variable
factor = None
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])