# -*- coding: utf-8 -*-
"""
Created on Tue Dec 23 17:16:57 2014

@author: jourdain
"""

import numpy as np
import pylab as plb
from scipy.stats import norm

#paramètres
T=1.#horizon
sig=0.2#volatilité
r=0.05#intérêt
S0=100.#valeur initiale du sous-jacent
K=110# Strike du Put considécontrré
N=1#nombre initial de pas de discrétisation
R=10#nombre de runs indépendants
M=1000000#taille des vecteurs de simulations indépendantes
RM=R*M#nombre total de simulations indépendantes
P=4#nombre d'itérations sur la valeur du pas




erreul=np.zeros(P)#vecteur des erreurs faibles Romberg Euler
liceul=np.zeros(P)#demi-largeur des intervalles de confiance 95% Romberg Euler
errmil=np.zeros(P)#vecteur des erreurs faibles Romberg Milstein
licmil=np.zeros(P)#demi-largeur des intervalles Romberg Milstein
Npas=np.zeros(P)# vecteur des nombres de pas



for i in range(P): #boucle sur le nombre de pas
  #Paramètres utiles pour les schémas avec N et 2N pas
  ######################################################
  ## À compléter avec le calcul des paramètres utiles
  #######################################################
  dt=T/N
  sigdt=sig*np.sqrt(dt)
  rdt=r*dt
  der=rdt-sigdt**2/2
  
  payeul=0
  careul=0
  paymil=0
  carmil=0
  for j in range(R):
      S=S0*np.ones(M)#initialisation Black-Scholes
      Se=S0*np.ones(M)#initialisation Eulen N pas
      S2e=S0*np.ones(M)#initialisation Eulen 2N pas
      Sm=S0*np.ones(M)#initialisation Milstein N pas
      S2m=S0*np.ones(M)#initialisation Milstein 2N pas
      for k in range(N): #boucle sur les pas de temps
          #génération d'un vecteur de gaussiennes centrées réduites comme somme de deux vecteurs
           g1=np.random.normal(0,1,M)/np.sqrt(2)
           g2=np.random.normal(0,1,M)/np.sqrt(2)
           g=g1+g2
     
          ############################################################
          ## À compléter avec l'evolution du sous-jacent S, 
          ## et des schémas Se,S2e et Sm,S2m
          ############################################################
           S=S*np.exp(sigdt*g+der)
           Se=Se*(1+sigdt*g+rdt)
           S2e=S2e*(1+sigdt*g1+rdt/2)*(1+sigdt*g2+rdt/2);
           Sm=Sm*(1+sigdt*g+(sigdt*g)**2/2+der)
           S2m=S2m*(1+sigdt*g1+(sigdt*g1)**2/2+der/2)*(1+sigdt*g2+(sigdt*g2)**2/2+der/2)

   #contribution des trajectoires des schémas d'Euler
      eul=2*(K-S2e)*(K>S2e)-(K-Se)*(K>Se)-(K-S)*(K>S)
      payeul=payeul+np.sum(eul)/M
      careul=careul+np.sum(eul**2)/M
      
      
  #contribution des trajectoires des schémas de Milstein
      mil=2*(K-S2m)*(K>S2m)-(K-Sm)*(K>Sm)-(K-S)*(K>S)
      paymil=paymil+np.sum(mil)/M
      carmil=carmil+np.sum(mil**2)/M
  erreul[i]=np.exp(-r*T)*payeul/R
  liceul[i]=1.96*np.exp(-r*T)*np.sqrt((careul/R-(payeul/R)**2)/RM)
  errmil[i]=np.exp(-r*T)*paymil/R
  licmil[i]=1.96*np.exp(-r*T)*np.sqrt((carmil/R-(paymil/R)**2)/RM)
  Npas[i]=N
  N=N*2#multiplication du nombre N de pas par 2

plb.plot(T/Npas**2,erreul, color="r", label="Romberg Euler")
plb.plot(T/Npas**2,erreul-liceul, color="b", label="IC 95%")
plb.plot(T/Npas**2,erreul+liceul, color="b")
plb.xlabel('Carre du pas de discretisation')
plb.legend(loc="best")
plb.show()

print("Romberg Euler")
print(erreul)
print("LIC Euler")
print(liceul)

plb.plot(T/Npas**2,errmil, color="r", label="Romberg 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()

print("Romberg Milstein")
print(errmil)
print("LIC Milstein")
print(licmil)