fig= plt.figure(figsize= (10, 5))
plt.plot(sigtested.index, data_return["Return"].values)
plt.show()

compare_hist_to_norm(data_return["Return"].values, bins= 25)

fig= plt.figure(figsize= (10, 5))
plt.plot(pickand_estimatored_k, c= "purple")
plt.title("Variation de l'estimateur de Pickands avec K")
plt.xlabel("k")
plt.ylabel("Pickland value")
plt.savefig('images/pickands_variation.png')
plt.show()

plt.figure(figsize= (10, 5))
plt.plot(rv_actual_r, "blue", label= "Log Return")
plt.plot(rv_predict_r, "red", label= "Prédiction")
plt.legend()
plt.title("Prédiction du log return à l'aide du modèle de Nuzman-Poor")
plt.savefig("images/predicted_actual_vol.png")
plt.show()

fig= plt.figure(figsize= (10, 5))
# Distribution de notre echantillon de backtest
plt.hist(backtest, bins= 30, color= "silver")

# VaR
plt.axvline(empiri_quant_esti, color= "red", label= "VaR empirique")
plt.axvline(parametric_quant_esti, color= "brown", label= "VaR paramétrique")
plt.axvline(nonparametric_quant_esti, color= "orange", label= "VaR non paramétrique")
plt.axvline(pickands_var_estimator_esti, color= "coral", label= "VaR Pickands")

# ES
plt.axvline(empiri_es_esti, color= "dodgerblue", label= "ES empirique")
plt.axvline(parametric_es_esti, color= "aqua", label= "ES paramétrique")
plt.axvline(nonparametric_es_esti, color= "skyblue", label= "ES non paramétrique")
plt.axvline(pickands_es_esti, color= "royalblue", label= "ES Pickands")

plt.legend(loc='best', frameon= True, shadow= True, framealpha= 0.9)
plt.title("Distribution de notre échantillon de test (avec {} observations)".format(len(backtest))\
          , size= 18)
plt.xlabel("Rendement")
plt.ylabel("Fréquence")

plt.savefig("images/var_es_backtest.png")

plt.show()

result_test

print("\t\t\t\t\t Distribution de la moyenne")
compare_hist_to_norm(mu_para, 50)
print("\n\t\t\t\t\t Distribution de l'écart-type")
compare_hist_to_norm(sig_para, 50)

print("\t\t\t\t\t Distribution de la VaR")
compare_hist_to_norm(var_para, 50)

print("\t\t\t\t\t Distribution de l'estimateur de Pickands")
compare_hist_to_norm(pick_, 30)
print("\n\t\t\t\t\t Distribution de l'extrème k")
compare_hist_to_norm(extreme_k_, 30)
print("\n\t\t\t\t\t Distribution de l'extrème k2")
compare_hist_to_norm(extreme_k2_, 30)
print("\n\t\t\t\t\t Distribution de l'extrème k4")
compare_hist_to_norm(extreme_k4_, 30)

print("\t\t\t\t\t Distribution de la VaR Pickands")
compare_hist_to_norm(var_pick_, 30)

obs_dist= stock_nonpara

fig = plt.figure(figsize=(12, 5))
ax = fig.add_subplot(111)

ax.hist(obs_dist, bins=20, density=True, label='Données',
        zorder=5, edgecolor='k', alpha=0.5)

for i in [0.0009, 0.0005, 0.0001, 0.00005]:
    kde = sm.nonparametric.KDEUnivariate(obs_dist)
    kde.fit(bw= i)
    ax.plot(kde.support, kde.density, lw=3, label='KDE avec h={}'.format(i), zorder=10)

plt.title("Distribution de notre VaR nonparamétrique")
ax.legend(loc='best')
ax.grid(True, zorder=-5)

for i in [0.99, 0.90, 0.95, 0.98, 0.995, 0.999, 0.9999]:
    val_min, idx_min = min((val, idx) for (idx, val) in enumerate(compute_all_es(i)))
    val_max, idx_max = max((val, idx) for (idx, val) in enumerate(compute_all_es(i)))
    print("Diamètre ES pour le seuil {}: {}".format(i, -val_max + val_min))
    
    val_min, idx_min = min((val, idx) for (idx, val) in enumerate(compute_all_var(i)))
    val_max, idx_max = max((val, idx) for (idx, val) in enumerate(compute_all_var(i)))
    print("Diamètre VaR pour le seuil {}: {}".format(i, -val_max + val_min))
    print("\n")

np.random.seed(1)

ampl= 1
data_noised= noised_dataset(amplitude= ampl)
fig= plt.figure(figsize= (10, 4))
plt.plot(data_noised["Return"])

noised_var= compute_all_var(seuil= 0.99, data= data_noised)
compute_diff_stab(amplitude= ampl, noised_dataset= noised_var)

np.random.seed(5)

ampl= 5
data_noised= noised_dataset(amplitude= ampl)
fig= plt.figure(figsize= (10, 4))
plt.plot(data_noised["Return"])

noised_var= compute_all_var(seuil= 0.99, data= data_noised)
compute_diff_stab(amplitude= ampl, noised_dataset= noised_var)

np.random.seed(10)

ampl= 10
data_noised= noised_dataset(amplitude= ampl)
fig= plt.figure(figsize= (10, 4))
plt.plot(data_noised["Return"])

noised_var= compute_all_var(seuil= 0.99, data= data_noised)
compute_diff_stab(amplitude= ampl, noised_dataset= noised_var)

np.random.seed(100)

ampl= 100
data_noised= noised_dataset(amplitude= ampl)
fig= plt.figure(figsize= (10, 4))
plt.plot(data_noised["Return"])

noised_var= compute_all_var(seuil= 0.99, data= data_noised)
compute_diff_stab(amplitude= ampl, noised_dataset= noised_var)

np.random.seed(100)

ampl= 100
data_noised= noised_dataset(amplitude= ampl)

w = pywt.Wavelet('haar')
x= data_noised["Return"].values

coeffs = pywt.wavedec(x, w, level= 1)
# EN enlevant le premier niveau du signal, on remarque que la diff de VaR est beaucoup plus failbe
# surtout pour la VaR historique
test= pywt.waverec(coeffs[:-1] + [None] * 1, w)#[:-6]
not_noised_data= pd.DataFrame(test, columns= ["Return"])

fig= plt.figure(figsize= (15, 5))

# print("\t\t\t\t\tPour le premier niveau de projection \n")
not_noised_var= compute_all_var(seuil= 0.99, data= not_noised_data)

plt.plot(x, label= "Noisy log return")
plt.plot(not_noised_data, label= "Denoised log return scale 1")
plt.title("Pour le premier niveau de projection")
plt.legend()
plt.show()
# compute_diff_stab(amplitude= ampl, noised_dataset= not_noised_var)

# sans bruit
fig= plt.figure(figsize= (15, 5))
plt.plot(data_return["Return"].values, color= "orange")
plt.title("Dataset initial (sans bruit)")
plt.show()

np.random.seed(100)

ampl= 100
data_noised= noised_dataset(amplitude= ampl)

w = pywt.Wavelet('haar')
x= data_noised["Return"].values

coeffs = pywt.wavedec(x, w, level= 3)
# EN enlevant le premier niveau du signal, on remarque que la diff de VaR est beaucoup plus failbe
# surtout pour la VaR historique
test1= pywt.waverec(coeffs[:-3] + [None] * 3, w)[:-6]
not_noised_data1= pd.DataFrame(test1, columns= ["Return"])

test2= pywt.waverec(coeffs[:-2] + [None] * 2, w)[:-6]
not_noised_data2= pd.DataFrame(test2, columns= ["Return"])

test3= pywt.waverec(coeffs[:-1] + [None] * 1, w)[:-6]
not_noised_data3= pd.DataFrame(test3, columns= ["Return"])

# notre signal bruité
# test4= pywt.waverec(coeffs[:4], w)[:-6]
# not_noised_data4= pd.DataFrame(test4, columns= ["Return"])

fig= plt.figure(figsize= (15, 5))

print("Pour le premier niveau de projection \n")
not_noised_var1= compute_all_var(seuil= 0.99, data= not_noised_data1)
compute_diff_stab(amplitude= ampl, noised_dataset= not_noised_var1)
print("\n\n")
print("Pour le deuxième niveau de projection \n")
not_noised_var2= compute_all_var(seuil= 0.99, data= not_noised_data2)
compute_diff_stab(amplitude= ampl, noised_dataset= not_noised_var2)
print("\n\n")
print("Pour le troisième niveau de projection \n")
not_noised_var3= compute_all_var(seuil= 0.99, data= not_noised_data3)
compute_diff_stab(amplitude= ampl, noised_dataset= not_noised_var3)

plt.plot(x, label= "Noisy log return")
plt.plot(not_noised_data3, label= "Denoised log return scale 3")
plt.plot(not_noised_data2, label= "Denoised log return scale 2")
plt.plot(not_noised_data1, label= "Denoised log return scale 1", color= "orange")
plt.title("Série bruité avec nos 3 niveaux de scale")
plt.legend()
plt.show()


# sans bruit
fig= plt.figure(figsize= (15, 5))
plt.plot(data_return["Return"].values, color= "skyblue")
plt.title("Dataset initial (sans bruit)")
plt.show()

%matplotlib inline

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl
import seaborn as sns
from sklearn.neighbors import KernelDensity
from sklearn.linear_model import LinearRegression
import random

import statsmodels.api as sm
from scipy import stats
import scipy.stats
import scipy.special as scsp
from scipy.stats import genextreme
import math

import pywt

import warnings
warnings.filterwarnings("ignore")

mpl.style.use("seaborn")

data= pd.read_csv("^GSPC_3.csv")
# data["Return"]= data["Close"].diff() / data["Close"]
data_copy= data.copy()
data["Return"]= data["Close"].shift(-1) / data["Close"]

data["Return"]= data["Return"].apply(np.log)
data.shape

# La gueule de notre série
plt.plot(data["Return"].values)
plt.title("Evolution du log rendement")
plt.savefig("images/data_return.png")
plt.show

data_return= data[["Date", "Return"]].dropna()

# Question 1

# Empirical quantile
empiri_quant= np.percentile(data_return["Return"], 1)
print("Empirical quantile : {}".format(empiri_quant))

# Paramétrique
alpha_01= scipy.stats.norm.ppf(0.01)
mu= np.mean(data_return["Return"], axis= 0)
# mu= np.mean(vol["Mean"])
sigma= np.std(data_return["Return"], axis= 0)
# sigma= np.mean(vol["sig"])
parametric_quant= alpha_01*sigma + mu
print("\nParametric quantile : {}".format(parametric_quant))

# Nonparametric distribution
# Monte Carlo simulation
np.random.seed(42) # pour fixer l'aléatoire car on fait une simulation
simulations = 1000000
sim_returns = np.random.normal(mu, sigma, simulations)
nonparametric_quant= np.percentile(sim_returns, 1)
# plt.hist(sim_returns, bins= 100)
print("\nNonparametric distribution : {}".format(nonparametric_quant))

# Question 2

# We define k = 10, to get the 5 worst returns

def pickands_estimator_func(data_return= data_return, k= 5):
    
    k= 5
    neg_data_return= - data_return["Return"].sort_values(ascending= False).values
    n= data_return.shape[0]

    extreme_return_k= max(neg_data_return[:n-k+1])
    extreme_return_k2= max(neg_data_return[:n-2*k+1])
    extreme_return_k4= max(neg_data_return[:n-4*k+1])



    return (1/np.log(2))*np.log((extreme_return_k - extreme_return_k2)\
                                         /(extreme_return_k2 - extreme_return_k4))


k= 5
neg_data_return= - data_return["Return"].sort_values(ascending= False).values
n= data_return.shape[0]

extreme_return_k= max(neg_data_return[:n-k+1])
extreme_return_k2= max(neg_data_return[:n-2*k+1])
extreme_return_k4= max(neg_data_return[:n-4*k+1])



pickands_estimator= (1/np.log(2))*np.log((extreme_return_k - extreme_return_k2)\
                                         /(extreme_return_k2 - extreme_return_k4))
pickands_estimator= pickands_estimator_func()
print("Pickands Estimator : {}".format(pickands_estimator))

pickand_estimatored_k= []

neg_data_return= - data_return["Return"].sort_values(ascending= False).values
n= data_return.shape[0]

for k in np.arange(1, 100):
    extreme_returned_k= max(neg_data_return[:n-k+1])
    extreme_returned_k2= max(neg_data_return[:n-2*k+1])
    extreme_returned_k4= max(neg_data_return[:n-4*k+1])



    pickand_estimatored_k.append((1/np.log(2))*np.log((extreme_returned_k - extreme_returned_k2)\
                                             /(extreme_returned_k2 - extreme_returned_k4)))
plt.plot(pickand_estimatored_k, c= "purple")
plt.title("Variation de l'estimateur de Pickands avec K")
plt.xlabel("k")
plt.ylabel("Pickland value")
plt.savefig('images/pickands_variation.png')
plt.show()

# Question 3

k= 5
p=0.99
coef= (pow((k/(n*(1-p))), pickands_estimator)-1)/(1-pow(2, -pickands_estimator))
pickands_var_estimator= -(coef*(extreme_return_k - extreme_return_k2) + extreme_return_k) # j'ai mis un - devant
                                                                                          # pour revenir au niveau d'une perte
print("La VaR 99%,grâce à la méthode des valeurs extrèmes est : {}"\
      .format(pickands_var_estimator))

# Question 4

b= 10
k= int(n/b) # nécessaire de round
seuil= np.log(0.03) # on définit notre seuil perte à 3%
#la valeur est bien positive car je multiplie mon dataset par -1

neg_value= - data_return["Return"].values #on garde l'ordre de notre dataset
somme_seuil= sum((neg_value > seuil))

# pour chaque partition (y en a k), je regarde si la valeur max dépasse mon seuil critique
somme_bucket_seuil= []
for i in np.arange(1, k+1):
    a= (i-1)*b
    z= i*b
    m= neg_value[a-1:z]
    if len(m) != 0:
        # dans chacun des buckets, je vérifie si le max dépasse mon seuil
        somme_bucket_seuil.append(max(neg_value[a-1:z]) > seuil)

extremal_index= sum(somme_bucket_seuil)/somme_seuil
print("Voici la valeur de l'indice extrêmal de Leadbetter : {}".format(extremal_index))

# La gueule de notre série
plt.plot(data["Return"].values)
plt.axhline(-0.03, color= "red")
plt.title("Evolution du log rendement")
plt.savefig("images/leadbetter_value.png")
plt.show()

sum(somme_bucket_seuil)

# Question 6 
# La seule méthode que nous avons trouvé qui permettait d'utiliser l'exposant de Hurst et la autocorrélation 
    # du mouvement brwonien fractionnaire est la prédiction de la log volatilité
    
# Calculons dans un premier temps un suite de log(return) mais pour plusieurs échelles de tps
data_copy= data.copy()

def return_per_i(i, data_copy= data_copy):
    data_copy["Return"+str(i)]= (data_copy["Close"].shift(-i) / data_copy["Close"]).apply(np.log)
    stocked= data_copy["Return"+str(i)].values[:-i].tolist()
    return stocked

# stock de notre valeur a régresser
final_stock= []
m_stocked= []
for i in [1, 2, 7, 30, 90, 100, 200, 252, 360]:
    
    stocked= return_per_i(i)
    n1= len(stocked)
    
    hurst_regression_data= []
    m_data= []
    for m in np.arange(1,len(stocked)):

        m_data.append(m)
        k= int(n1/m)
        agg_return= []
        for i in np.arange(1, k+1):
            t= 1+m*(i-1)
            agg_return.append(np.sum(stocked[t:k*m])/m)

        hurst_regression_data.append(np.sum(np.abs(agg_return[0: k]))/k)
        
    final_stock += hurst_regression_data
    m_stocked += m_data

# x_log= np.log(np.arange(1, len(stocked)))
# x_log= np.log(np.arange(1, 668).tolist()*6)[:-3]
x_log= np.log(m_stocked)
intercept= np.repeat(1, len(m_stocked))
x= np.column_stack((intercept, x_log))
y= final_stock
lr= LinearRegression(fit_intercept= True).fit(x, y)

h= lr.coef_[1] + 1 #mon exposant de Hurst
h

sigtested= data[["Date", "Close"]]
sigtested["Close"]= np.log(data["Close"])
sigtested["Date"]= pd.to_datetime(sigtested["Date"], format= "%Y-%m-%d")
sigtested.set_index(["Date"], inplace= True)

sigtested["Return"]= sigtested["Close"].shift(-1) / sigtested["Close"]
sigtested= sigtested[["Return"]].dropna()

def c_tilde(h):
    return scsp.gamma(3. / 2. - h) / scsp.gamma(h + 1. / 2.) * scsp.gamma(2. - 2. * h)

def forecast_XTS(rvdata, h, date, nLags, delta, nu):
    i = np.arange(nLags)
    cf = 1./((i + 1. / 2.) ** (h + 1. / 2.) * (i + 1. / 2. + delta))
    ldata = rvdata.truncate(after=date)
    # ldata= rvdata
    l = len(ldata)
    ldata = np.log(ldata.iloc[l - nLags:])
    ldata['cf'] = np.fliplr([cf])[0]

    ldata = ldata.dropna()
    fcst = (ldata.iloc[:, 0] * ldata['cf']).sum() / sum(ldata['cf'])

    return math.exp(fcst + 2 * nu ** 2 * c_tilde(h) * delta ** (2 * h))

x = np.arange(1, 100)

def dlsig2(sic, x, pr=False):
    if pr:
        a= np.array([(sig-sig.shift(lag)).dropna() for lag in x])
        a=a ** 2
    return [np.mean((sig-sig.shift(lag)).dropna() ** 2) for lag in x]


nu_= []
sig = sigtested
sig =  np.log(np.sqrt(sig))
sig = sig.dropna()
model = np.polyfit(np.log(x), np.log(dlsig2(sig, x)), 1)
nu_.append(np.sqrt(np.exp(model[1])))

nu= nu_[0]
rvdata= sigtested

n = len(rvdata)
delta = 1
nLags = 300
dates = rvdata.iloc[nLags:n-delta].index
rv_predict = [forecast_XTS(rvdata, h=h, date=d, nLags=nLags,
                           delta=delta, nu=nu) for d in dates]
rv_actual = rvdata.iloc[nLags+delta:n].values

nu= nu_[0]
rvdata= sigtested

n = len(rvdata)
delta = 1
nLags = 300
dates = rvdata.iloc[nLags:n-delta].index
rv_predict = [forecast_XTS(rvdata, h=h, date=d, nLags=nLags,
                           delta=delta, nu=nu) for d in dates]
rv_actual = rvdata.iloc[nLags+delta:n].values

dates
new_date= ['2020-02-12', '2020-02-13', '2020-02-14', '2020-02-15', '2020-02-16'
           , '2020-02-17', '2020-02-18', '2020-02-19', '2020-02-20', '2020-02-21'
           , '2020-02-22', '2020-02-23', '2020-02-24', '2020-02-25', '2020-02-26']

rv_predict_test = [forecast_XTS(rvdata, h=h, date=d, nLags=nLags,
                           delta=delta, nu=nu) for d in new_date]

# Prédiction du return en rouge
rv_actual_r = (rvdata.iloc[nLags+delta:n].shift(-1) - rvdata.iloc[nLags+delta:n]).values
rv_predict_r= pd.DataFrame(rv_predict, columns= ["Return"])
rv_predict_r= (rv_predict_r["Return"].shift(-1) - rv_predict_r["Return"]).values

plt.plot(rv_actual_r, "blue")
plt.plot(rv_predict_r, "red")
plt.savefig("images/predicted_actual_vol.png")
plt.show()

# Question 7
# A partir d'une VaR on peut calculer une ES. 
# Un des problèmes de la VaR c'est qu'on a aucune idée de ce qu'il
# se passe outre passé notre dite-VaR. L'expected shortfall est la moyenne
# des rendements une fois passé notre VaR. Ainsi, on garde l'information 
# de la gueule de notre queue après la VaR.

# On va donc calculer la moyenne aprèsn notre VaR

# ES empirique
empiri_quant
empiri_es= np.mean(data_return[data_return["Return"] < empiri_quant]["Return"])
print(empiri_es)

# # ES parametrique
parametric_quant
parametric_es= np.mean(data_return[data_return["Return"] < parametric_quant]["Return"])
print(parametric_es)


# # ES non-parametrique
nonparametric_quant
nonparametric_es= np.mean(data_return[data_return["Return"] < nonparametric_quant]["Return"])
print(nonparametric_es)

# # ES Pickands
pickands_var_estimator
pickands_es= np.mean(data_return[data_return["Return"] < pickands_var_estimator]["Return"])
print(pickands_es)

# Question 8

from sklearn.model_selection import train_test_split
np.random.seed(55)
estimation, backtest= train_test_split(data_return["Return"], test_size= 0.3)

# Empirical quantile
empiri_quant_esti= np.percentile(estimation, 1)

# Parametric quantile
# alpha_01= -3.090232
# vol_esti= pd.DataFrame()
# vol_esti["Return"]= estimation.copy()
# vol_esti["Return_t1"]= vol_esti["Return"].shift(-1)
# vol_esti["Mean"]= (vol_esti["Return"] + vol_esti["Return_t1"])/2
# vol_esti["sig"]= np.sqrt(pow(vol_esti["Return"] - vol_esti["Mean"], 2) + pow(vol_esti["Return_t1"] - vol_esti["Mean"], 2))

alpha_01= scipy.stats.norm.ppf(0.01)
mu= np.mean(estimation, axis= 0)
# mu_esti= np.mean(vol_esti["Mean"])
sigma= np.std(estimation, axis= 0)
# sigma_esti= np.mean(vol_esti["sig"])
parametric_quant_esti= alpha_01*sigma + mu

# Nonparametric distribution
np.random.seed(42)
simulations = 1000000
sim_returns = np.random.normal(mu, sigma, simulations)
nonparametric_quant_esti= np.percentile(sim_returns, 1)

# Pickands VaR
k= 5
neg_data_return= - estimation.sort_values(ascending= False).values
n= estimation.shape[0]

extreme_return_k= max(neg_data_return[:n-k+1])
extreme_return_k2= max(neg_data_return[:n-2*k+1])
extreme_return_k4= max(neg_data_return[:n-4*k+1])

pickands_estimator_esti= (1/np.log(2))*np.log((extreme_return_k - extreme_return_k2)\
                                         /(extreme_return_k2 - extreme_return_k4))


p=0.99
coef= (pow((k/(n*(1-p))), pickands_estimator_esti)-1)/(1-pow(2, -pickands_estimator_esti))
pickands_var_estimator_esti= -(coef*(extreme_return_k - extreme_return_k2) + extreme_return_k) # j'ai mis un - devant
                                                                                          # pour revenir au niveau d'une perte

# ES

# ES empirique
empiri_es_esti= np.mean(estimation[estimation < empiri_quant_esti])

# # ES parametrique
parametric_es_esti= np.mean(estimation[estimation < parametric_quant_esti])


# # ES non-parametrique
nonparametric_es_esti= np.mean(estimation[estimation < nonparametric_quant_esti])


# # ES Pickands
pickands_es_esti= np.mean(estimation[estimation < pickands_var_estimator_esti])

fig= plt.figure(figsize=(10, 7))

stock_value= {
    "VaR empirique":empiri_quant_esti,
    "VaR paramétrique":parametric_quant_esti,
    "VaR Non-paramétrique":nonparametric_quant_esti,
    "VaR Pickands":pickands_var_estimator_esti,
    "ES empirique":empiri_es_esti,
    "ES paramétrique":parametric_es_esti,
    "ES Non-paramétrique":nonparametric_es_esti,
    "ES Pickands":pickands_es_esti,
}

# Distribution de notre echantillon de backtest
plt.hist(backtest, bins= 30, color= "silver")

# VaR
plt.axvline(empiri_quant_esti, color= "red", label= "VaR empirique")
plt.axvline(parametric_quant_esti, color= "brown", label= "VaR paramétrique")
plt.axvline(nonparametric_quant_esti, color= "orange", label= "VaR non paramétrique")
plt.axvline(pickands_var_estimator_esti, color= "coral", label= "VaR Pickands")

# ES
plt.axvline(empiri_es_esti, color= "dodgerblue", label= "ES empirique")
plt.axvline(parametric_es_esti, color= "aqua", label= "ES paramétrique")
plt.axvline(nonparametric_es_esti, color= "skyblue", label= "ES non paramétrique")
plt.axvline(pickands_es_esti, color= "royalblue", label= "ES Pickands")

plt.legend(loc='best', frameon= True, shadow= True, framealpha= 0.9)
plt.title("Distribution de notre échantillon de test (avec {} observations)".format(len(backtest))\
          , size= 18)
plt.xlabel("Rendement")
plt.ylabel("Fréquence")

plt.savefig("images/var_es_backtest.png")

plt.show()

# backtesting : Conditional Test by Acerbi and Szekely

num_hit= sum(backtest < nonparametric_quant_esti) #nombre de fois ou mon rendement est en dessous de mon seuil
# if num_hist > 0:
    
# else:
#     print("Pas possible")
# np.sum(((backtest[backtest < nonparametric_quant_esti] / nonparametric_es_esti) + 1)) / num_hit

def couverture_nonconditionnelle(seuil, data= backtest):
    """
    C'est la proba qu'un rendement soit inférieur à mon seuil
    """
    return sum(data < seuil)/len(data)

def test_couverture_nonconditionnelle(seuil, data= backtest, p= 0.01):
    """
    On cherche à savoir si la proba de couverture non conditionnelle est statistiquement
        différente de 0. On fait donc que le ratio num_hit/n tend vers notre proba de risque alpha
        et que le nombre total de num_hit suit une loi binomiale.
        
        On va comparer notre valeur compute au quantile de la loi normale centrée réduite
    """
    num_hit= sum(data < seuil)
    n= len(data)
    return (num_hit - p*n)/np.sqrt(p*(1-p)*n) # suit suivre une loi normale centrée réduite


def test_normality_nonconditionnelle(seuil, data= backtest, p= 0.01):
    import scipy.stats
    
    if np.abs(test_couverture_nonconditionnelle(seuil= seuil, data= data, p= p)) < scipy.stats.norm.ppf(1-p):
        # notre mesure du risque est "bonne" (on sur-estime le risque)
#         print("On accepte l'hypothèse H0 : notre probabilité de hit est égale à notre seuil alpha= {}".format(p))
        return "On accepte H0"
    else:
#         print("On rejète l'hypothèse H0")
        return "On rejète H0"

# couverture_nonconditionnelle(nonparametric_quant_esti)
# test_couverture_nonconditionnelle(nonparametric_quant_esti)
test_normality_nonconditionnelle(nonparametric_quant_esti)
test_normality_nonconditionnelle(pickands_es_esti)

def stat_kupiec(seuil, data= backtest, p= 0.01):
    num_hit= sum(data < seuil)
    n= len(data)
    minus= -2*np.log(pow(1-p, n - num_hit)*pow(p, num_hit))
    plus= 2*np.log(pow(1- (num_hit/n), n - num_hit)*pow(num_hit, num_hit))
    return minus + plus

def test_chi_nonconditionnelle(seuil, data= backtest, p= 0.01, k= 1): # 1 degree de liberté
    import scipy.stats
    
    if stat_kupiec(seuil= seuil, data= data, p= p) < scipy.stats.chi2.ppf(1-p, df= k):
        # notre mesure du risque est "bonne" (on sur-estime le risque)
#         print("On accepte l'hypothèse H0 : notre probabilité de hit est égale à notre seuil alpha= {}".format(p))
        return "On accepte H0"
    else:
#         print("On rejète l'hypothèse H0")
        return "On rejète H0"

stock_chi= {
    "VaR empirique":test_chi_nonconditionnelle(empiri_quant_esti),
    "VaR paramétrique":test_chi_nonconditionnelle(parametric_quant_esti),
    "VaR Non-paramétrique":test_chi_nonconditionnelle(nonparametric_quant_esti),
    "VaR Pickands":test_chi_nonconditionnelle(pickands_var_estimator_esti),
    "ES empirique":test_chi_nonconditionnelle(empiri_es_esti),
    "ES paramétrique":test_chi_nonconditionnelle(parametric_es_esti),
    "ES Non-paramétrique":test_chi_nonconditionnelle(nonparametric_es_esti),
    "ES Pickands":test_chi_nonconditionnelle(pickands_es_esti),
}
stock_norm= {
    "VaR empirique":test_normality_nonconditionnelle(empiri_quant_esti),
    "VaR paramétrique":test_normality_nonconditionnelle(parametric_quant_esti),
    "VaR Non-paramétrique":test_normality_nonconditionnelle(nonparametric_quant_esti),
    "VaR Pickands":test_normality_nonconditionnelle(pickands_var_estimator_esti),
    "ES empirique":test_normality_nonconditionnelle(empiri_es_esti),
    "ES paramétrique":test_normality_nonconditionnelle(parametric_es_esti),
    "ES Non-paramétrique":test_normality_nonconditionnelle(nonparametric_es_esti),
    "ES Pickands":test_normality_nonconditionnelle(pickands_es_esti),
}

result_test= pd.DataFrame()
result_test["Méthode"]= stock_chi.keys()
result_test["Valeur"]= result_test["Méthode"].map(stock_value)
result_test["Test non conditionnel"]= result_test["Méthode"].map(stock_norm)
result_test["Test Kupiec"]= result_test["Méthode"].map(stock_chi)

print(result_test)

# Question 9

# Question 10

# On va faire N tirage de  avec remise dans notre dataset de 10 log return, puis en faire la moyenne
# et l'écart type
random.seed(55)

mu_para= []
sig_para= []
var_para= []

all_return= data_return["Return"].values
k= 10
alpha_01= scipy.stats.norm.ppf(0.01)

for i in range(100000):
    random_values= random.choices(all_return, k= k)
    
    mu_test= np.mean(random_values)
    sig_test= np.std(random_values)
    
    var_para.append(alpha_01*sig_test + mu_test)
    mu_para.append(mu_test)
    sig_para.append(sig_test)
    

def compare_hist_to_norm(data, bins= 25):
    fig= plt.figure(figsize= (10, 5))
    from scipy.stats import skewnorm

#     mu, std, sk, k= scipy.stats.norm.stats(data, moments='mvsk')
    mu, std = scipy.stats.norm.fit(data)
#     params= scipy.stats.norm.fit(data)
#     sk= scipy.stats.skewnorm(data)
    # Plot l'hist
    plt.hist(data, bins= bins, density=True, alpha=0.6, color='purple', label= "Données")

    # Plot le PDF.
    xmin, xmax = plt.xlim()
    X = np.linspace(xmin, xmax)
    
    plt.plot(X, scipy.stats.norm.pdf(X, mu, std), label= "Normal Distribution")
    plt.plot(X, skewnorm.pdf(X, *skewnorm.fit(data)), color= 'black', label= "Skewed Normal Distribution")
    
    mu, std = scipy.stats.norm.fit(data)
    sk= scipy.stats.skew(data)
    
    title2 = "Moments mu: {}, sig: {}, sk: {}".format(round(mu, 4), round(std, 4), round(sk, 4))
    plt.ylabel("Fréquence", rotation= 90)
    plt.title(title2)
    plt.legend()
    plt.show()

compare_hist_to_norm(mu_para, 50)
compare_hist_to_norm(sig_para, 50)
compare_hist_to_norm(var_para, 50)

# Question 11
# si le prix suit une gaussienne iid, alors on va générer plusieurs distribution de prix en
# prenant compte de la moyenne et de la volatilité de mes prix
# Attention, nous travaillons sur des log(prix)

data_second= data_return["Return"]
# les 2 premiers moments de mes log prix
mean_lp= np.mean(data_second, axis= 0)
sig_lp= np.std(data_second, axis= 0)

pick_= []
extreme_k_= []
extreme_k2_= []
extreme_k4_= []
var_pick_= []
for i in range(10000):
    simulations= 1000 # on va générer 1000 prix gaussian iid
    log_prix= pd.DataFrame(np.random.normal(mean_lp, sig_lp, simulations), columns= ["Close"])
    log_rendement= pd.DataFrame(np.random.normal(mean_lp, sig_lp, simulations), columns= ["Return"])
    
    # calcul var evt
    returns_= pd.DataFrame(log_rendement["Return"].values, columns= ["Return"]).dropna()
    
    k= 5
    neg_data_return= - returns_["Return"].sort_values(ascending= False).values
    n= returns_.shape[0]
    extreme_return_k= max(neg_data_return[:n-k+1])
    extreme_return_k2= max(neg_data_return[:n-2*k+1])
    extreme_return_k4= max(neg_data_return[:n-4*k+1])

    pickands_estimator= (1/np.log(2))*np.log((extreme_return_k - extreme_return_k2)\
                                             /(extreme_return_k2 - extreme_return_k4))
    
    
    k= 5
    p=0.99
    coef= (pow((k/(n*(1-p))), pickands_estimator)-1)/(1-pow(2, -pickands_estimator))
    pickands_var_estimator= -(coef*(extreme_return_k - extreme_return_k2) + extreme_return_k)
    
    
    pick_.append(pickands_estimator)
    extreme_k_.append(extreme_return_k)
    extreme_k2_.append(extreme_return_k2)
    extreme_k4_.append(extreme_return_k4)

    var_pick_.append(pickands_var_estimator)

compare_hist_to_norm(pick_, 30)

compare_hist_to_norm(extreme_k_, 30)

compare_hist_to_norm(extreme_k2_, 30)

compare_hist_to_norm(extreme_k4_, 30)

compare_hist_to_norm(var_pick_, 30)

# Question 12

# Question 13

mu= np.mean(data_return["Return"].values)
sigma= np.std(data_return["Return"].values)

stock_nonpara= []
for i in range(1000):
    simulations = 1000
    sim_returns = np.random.normal(mu, sigma, simulations)
    nonparametric_quant= np.percentile(sim_returns, 1)
    stock_nonpara.append(nonparametric_quant)

obs_dist= stock_nonpara

fig = plt.figure(figsize=(12, 5))
ax = fig.add_subplot(111)

ax.hist(obs_dist, bins=20, density=True, label='Données',
        zorder=5, edgecolor='k', alpha=0.5)

for i in [0.0009, 0.0005, 0.0001, 0.00005]:
    kde = sm.nonparametric.KDEUnivariate(obs_dist)
    kde.fit(bw= i)
    ax.plot(kde.support, kde.density, lw=3, label='KDE avec h={}'.format(i), zorder=10)


ax.legend(loc='best')
ax.grid(True, zorder=-5)

# QUestion 14

def compute_all_var(seuil= 0.99, data= data_return):
    
    estimation= data["Return"]

    # Empirical quantile
    empiri_quant_esti= np.percentile(estimation, (1-seuil)*100)

    # Parametric quantile
    # alpha_01= -3.090232
    # vol_esti= pd.DataFrame()
    # vol_esti["Return"]= estimation.copy()
    # vol_esti["Return_t1"]= vol_esti["Return"].shift(-1)
    # vol_esti["Mean"]= (vol_esti["Return"] + vol_esti["Return_t1"])/2
    # vol_esti["sig"]= np.sqrt(pow(vol_esti["Return"] - vol_esti["Mean"], 2) + pow(vol_esti["Return_t1"] - vol_esti["Mean"], 2))    
    alpha_01= scipy.stats.norm.ppf(1-seuil)
    mu= np.mean(estimation, axis= 0)
    # mu_esti= np.mean(vol_esti["Mean"])
    sigma= np.std(estimation, axis= 0)
    # sigma_esti= np.mean(vol_esti["sig"])
    parametric_quant_esti= alpha_01*sigma + mu

    # Nonparametric distribution
#     np.random.seed(42)
    simulations = 1000000
    sim_returns = np.random.normal(mu, sigma, simulations)
    nonparametric_quant_esti= np.percentile(sim_returns, (1-seuil)*100)

    # Pickands VaR
    k= 5
    neg_data_return= - estimation.sort_values(ascending= False).values
    n= estimation.shape[0]

    extreme_return_k= max(neg_data_return[:n-k+1])
    extreme_return_k2= max(neg_data_return[:n-2*k+1])
    extreme_return_k4= max(neg_data_return[:n-4*k+1])

    pickands_estimator_esti= (1/np.log(2))*np.log((extreme_return_k - extreme_return_k2)\
                                             /(extreme_return_k2 - extreme_return_k4))


    p= seuil
    coef= (pow((k/(n*(1-p))), pickands_estimator_esti)-1)/(1-pow(2, -pickands_estimator_esti))
    pickands_var_estimator_esti= -(coef*(extreme_return_k - extreme_return_k2) + extreme_return_k) # j'ai mis un - devant
                                                                                              # pour revenir au niveau d'une perte
    return (empiri_quant_esti, parametric_quant_esti, nonparametric_quant_esti, pickands_var_estimator_esti)

def compute_all_es(seuil= 0.99):
    
    estimation= data_return["Return"]
    res= compute_all_var(seuil)
    # ES empirique
    empiri_es_esti= np.mean(estimation[estimation < res[0]])

    # # ES parametrique
    parametric_es_esti= np.mean(estimation[estimation < res[1]])

    # # ES non-parametrique
    nonparametric_es_esti= np.mean(estimation[estimation < res[2]])

    # # ES Pickands
    pickands_es_esti= np.mean(estimation[estimation < res[3]])
    
    return [empiri_es_esti, parametric_es_esti, nonparametric_es_esti, pickands_es_esti]

i= 0.99
val_min, idx_min = min((val, idx) for (idx, val) in enumerate(compute_all_es(i)))
val_max, idx_max = max((val, idx) for (idx, val) in enumerate(compute_all_es(i)))
print("Diameter ES pour le seuil {}: {}".format(i, -val_max + val_min))

val_min, idx_min = min((val, idx) for (idx, val) in enumerate(compute_all_var(i)))
val_max, idx_max = max((val, idx) for (idx, val) in enumerate(compute_all_var(i)))
print("Diameter VaR pour le seuil {}: {}".format(i, -val_max + val_min))
print("\n")

# Question 15
for i in [0.99, 0.90, 0.95, 0.98, 0.995, 0.999, 0.9999]:
    val_min, idx_min = min((val, idx) for (idx, val) in enumerate(compute_all_es(i)))
    val_max, idx_max = max((val, idx) for (idx, val) in enumerate(compute_all_es(i)))
    print("Diameter ES pour le seuil {}: {}".format(i, -val_max + val_min))
    
    val_min, idx_min = min((val, idx) for (idx, val) in enumerate(compute_all_var(i)))
    val_max, idx_max = max((val, idx) for (idx, val) in enumerate(compute_all_var(i)))
    print("Diameter VaR pour le seuil {}: {}".format(i, -val_max + val_min))
    print("\n")
    
# certaines valeurs ne peuvant pas être calculés (ES) car pour un seuil donné, nous n'avons pas
# suffisamment de données hist.

# Dans tous les cas, l'écart est plus faible avec l'ES que la VaR ==> moins de variance entre
# les modèles ES que VaR

# Question 16
# np.random.seed(55)
def noised_dataset(amplitude= 1):
    val= data["Close"]

    noise= pd.DataFrame(np.random.normal(0, amplitude, len(val)).tolist(), columns= ["Noise"])
    data_noised= pd.DataFrame(val.values, columns= ["Close"])
    data_noised["Close_noised"]= data_noised["Close"] + noise["Noise"]
    # data_noised.head()
    data_noised["Return"]= np.log((data_noised["Close_noised"].shift(-1) / data_noised["Close_noised"]))
    data_noised= data_noised.dropna()
    
    return data_noised

data_noised= noised_dataset()
plt.plot(data_noised["Return"])
# archi bruité

good_var= compute_all_var()
noised_var= compute_all_var(seuil= 0.99, data= data_noised)

def compute_diff_stab(amplitude, noised_dataset, gv= good_var):
    type_var= ["Historique", "Paramétrique", "Non Paramétrique", "Pickands"]
    for i in np.arange(0, 4):
        diff= np.abs(np.abs(good_var[i]) - np.abs(noised_dataset[i]))
        print("Pour un bruit blanc gaussien avec une amplitude de {},\n la différence de valeur de la VaR {} est: {}\n"\
             .format(amplitude, type_var[i], diff))
    pass

ampl= 1
data_noised= noised_dataset(amplitude= ampl)
plt.plot(data_noised["Return"])

noised_var= compute_all_var(seuil= 0.99, data= data_noised)
compute_diff_stab(amplitude= ampl, noised_dataset= noised_var)
# pour un bruit blanc d'amplitude 1, c'est la VaR Pickands qui est la plus robuste

ampl= 5
data_noised= noised_dataset(amplitude= ampl)
plt.plot(data_noised["Return"])

noised_var= compute_all_var(seuil= 0.99, data= data_noised)
compute_diff_stab(amplitude= ampl, noised_dataset= noised_var)
# pour un bruit blanc d'amplitude 5, c'est la VaR Paramétrique qui est la plus robuste

ampl= 20
data_noised= noised_dataset(amplitude= ampl)
plt.plot(data_noised["Return"])

noised_var= compute_all_var(seuil= 0.99, data= data_noised)
compute_diff_stab(amplitude= ampl, noised_dataset= noised_var)
# pour un bruit blanc d'amplitude 20, c'est la VaR Historique qui est la plus robuste

ampl= 100
data_noised= noised_dataset(amplitude= ampl)
plt.plot(data_noised["Return"])

noised_var= compute_all_var(seuil= 0.99, data= data_noised)
compute_diff_stab(amplitude= ampl, noised_dataset= noised_var)
# pour un bruit blanc d'amplitude 100, c'est la VaR Paramétrique qui est la plus robuste

# Question 17

import pywt

ampl= 100
data_noised= noised_dataset(amplitude= ampl)

w = pywt.Wavelet('haar')
x= data_noised["Return"].values
# plt.plot(w.dec_lo)
coeffs = pywt.wavedec(x, w, level= 1)

def reconstruction_plot(yyy, **kwargs):
    """Plot signal vector on x [0,1] independently of amount of values it contains."""
    plt.plot(np.linspace(0, 1, len(yyy)), yyy, **kwargs)

def haarFWT ( signal, level ):

    s = .5;                  # scaling -- try 1 or ( .5 ** .5 )

    h = [ 1,  1 ];           # lowpass filter
    g = [ 1, -1 ];           # highpass filter        
    f = len ( h );           # length of the filter

    t = signal;              # 'workspace' array
    l = len ( t );           # length of the current signal
    y = [0] * l;             # initialise output

#     t = t + [ 0, 0 ];        # padding for the workspace

    for i in range ( level ):

        y [ 0:l ] = [0] * l; # initialise the next level 
        l2 = l // 2;         # half approximation, half detail

        for j in range ( l2 ):            
            for k in range ( f ):                
                y [j]    += t [ 2*j + k ] * h [ k ] * s;
                y [j+l2] += t [ 2*j + k ] * g [ k ] * s;

        l = l2;              # continue with the approximation
        t [ 0:l ] = y [ 0:l ] ;

    return y

plt.figure(figsize= (10, 5))
plt.plot(x)
plt.plot(haarFWT(x.tolist(), 1))

# plt.plot(data_return["Return"].values)
# plt.plot(haarFWT(x.tolist(), 2))
# plt.plot(haarFWT(x.tolist(), 3))
# plt.plot(haarFWT(x.tolist(), 4))
# plt.plot(haarFWT(x.tolist(), 5))
plt.legend(["Données", "notre lvl 1"])
# j'arrive pas a le développer moi même car ça ne marche pas

# EN enlevant le premier niveau du signal, on remarque que la diff de VaR est beaucoup plus failbe
# surtout pour la VaR historique
# test= data_noised["Return"].values - pywt.waverec(coeffs[:-3] + [None] * 3, w)[:-6]
test= pywt.waverec(coeffs[:-1] + [None] * 1, w)#[:-6]
not_noised_data= pd.DataFrame(test, columns= ["Return"])
plt.hist(test, bins= 30)
plt.show()

fig= plt.figure(figsize= (10, 5))
ampl= 10
print("Pour le premier niveau de projection \n")
noised_var= compute_all_var(seuil= 0.99, data= not_noised_data)
plt.plot(not_noised_data)
compute_diff_stab(amplitude= ampl, noised_dataset= noised_var)

# test= data_noised["Return"].values -  pywt.waverec(coeffs[:-2] + [None] * 2, w)[:-6]
# not_noised_data= pd.DataFrame(test, columns= ["Return"])

# ampl= 10
# print("Pour le deuxième niveau de projection \n")
# noised_var= compute_all_var(seuil= 0.99, data= not_noised_data)
# # plt.plot(not_noised_data)
# compute_diff_stab(amplitude= ampl, noised_dataset= noised_var)

# test= data_noised["Return"].values - pywt.waverec(coeffs[:-1] + [None] * 1, w)[:-2] 
# not_noised_data= pd.DataFrame(test, columns= ["Return"])

# ampl= 10
# print("Pour le troisième niveau de projection \n")
# noised_var= compute_all_var(seuil= 0.99, data= not_noised_data)
# # plt.plot(not_noised_data)
# compute_diff_stab(amplitude= ampl, noised_dataset= noised_var)

# ampl= 10
# data_noised= noised_dataset(amplitude= ampl)
# plt.plot(data_noised["Return"])

# noised_var= compute_all_var(seuil= 0.99, data= data_noised)
# compute_diff_stab(amplitude= ampl, noised_dataset= noised_var)