# -*- coding: utf-8 -*-
"""
Created on Tue Dec 23 16:24:10 2014

@author: jourdain
"""


import numpy as np
import pylab as plb


T=1. #horizon
sig=0.2 #volatilité
r=0.05 #taux d'intérêt
S0=100 #valeur initial du sous-jacent
N=1 #nombre initial de pas de discrétisation
M=100000 #nombre de simulations indépendantes
P=6#nombre d'itérations sur la valeur d pas

erreul=np.zeros(P)
liceul=np.zeros(P)
errmil=np.zeros(P)
licmil=np.zeros(P)
Npas=np.zeros(P)

for i in range(P):
    #Paramètres utiles pour la discrétisation avec N pas
    ##################################################
    ## À compléter avec le calcul des paramètres
    ##################################################
   
    
    S=S0*np.ones(M) #initialisation des M traj de l'EDS
    Se=S0*np.ones(M) #initialisation des M traj du schéma d'Euler
    Sm=S0*np.ones(M) #initialisation des M traj du schéma de Milstein
    maxerreul=np.zeros(M) #initialisation de l'écart maximum EDS-Euler en valeur absolue
    maxerrmil=np.zeros(M) #initialisation de l'écart maximum EDS-Milstein en valeur absolue
    for k in range(N):
        g=np.random.normal(0,1,M)#génération d'un vecteur de gaussiennes centrées réduites
        ############################################################
        ## À compléter avec l'evolution du sous-jacent S, du schéma d'Euler Se
        ## et du schéma de Milstein Sm
        ############################################################
        
        
        
        maxerreul=np.max([maxerreul,np.abs(S-Se)],axis=0)
        maxerrmil=np.max([maxerrmil,np.abs(S-Sm)],axis=0)
    someul=np.sum(maxerreul**2)
    careul=np.sum(maxerreul**4)
    erreul[i]=someul/M
    liceul[i]=1.96*np.sqrt((careul/M-(someul/M)**2)/M)
    sommil=np.sum(maxerrmil**2)
    carmil=np.sum(maxerrmil**4)
    errmil[i]=sommil/M
    licmil[i]=1.96*np.sqrt((carmil/M-(sommil/M)**2)/M)
    Npas[i]=N
    N=2*N  
plb.plot(T/Npas,erreul, color="r", label="Erreur quadratique Euler")
plb.plot(T/Npas,erreul-liceul, color="b", label="IC 95%")
plb.plot(T/Npas,erreul+liceul, color="b")
plb.xlabel('Pas de discretisation')
plb.legend(loc="best")
plb.show()
#Affichage des vecteurs d'erreurs et des demi-largeurs des IC à 95%
print("Erreur euler")
print(erreul)
print("LIC Euler")
print(liceul)
plb.plot(T/Npas**2,errmil, color="r", label="Erreur quadratique Milstein")
plb.plot(T/Npas**2,errmil-licmil, color="b", label="IC 95%")
plb.plot(T/Npas**2,errmil+licmil, color="b")
plb.xlabel('Carre du pas de discretisation')
plb.legend(loc="best")
plb.show()
#Affichage des vecteurs d'erreurs et des demi-largeurs des IC à 95%
print("Erreur Milstein")
print(errmil)
print("LIC Milstein")
print(licmil)