# -*- 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
M=100000#nombre de simulations indépendantes
P=6#nombre d'itérations sur la valeur du pas

#formule Black-Scholes pour le prix du Put
d1=(np.log(S0/K)+r*T)/(sig*np.sqrt(T))+sig*np.sqrt(T)/2
d2=d1-sig*np.sqrt(T)
BS=K*np.exp(-r*T)*norm.cdf(-d2)-S0*norm.cdf(-d1)


erreul=np.zeros(P)#vecteur des erreurs faibles Euler
liceul=np.zeros(P)#demi-largeur des intervalles de confiance 95% Euler
conterreul=np.zeros(P)#vecteur des erreurs faibles Euler avec variable de contrôle
contliceul=np.zeros(P)#demi-largeur des intervalles de confiance 95%
errmil=np.zeros(P)#vecteur des erreurs faibles Euler
licmil=np.zeros(P)#demi-largeur des intervalles de confiance 95% Euler
conterrmil=np.zeros(P)#vecteur des erreurs faibles Euler avec variable de contrôle
contlicmil=np.zeros(P)#demi-largeur des intervalles de confiance 95%
Npas=np.zeros(P)# vecteur des nombres de pas



for i in range(P): #boucle sur le nombre de pas
  #Paramètres utiles pour la discrétisation avec N pas
  dt=T/N
  sigdt=sig*np.sqrt(dt)
  rdt=r*dt
  der=rdt-sigdt**2/2
  
  S=S0*np.ones(M);
  Se=S0*np.ones(M);
  Sm=S0*np.ones(M);
  for k in range(N): #boucle sur les pas de temps
      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, 
      ## et des schémas Se et Sm
      ############################################################
      S=S*np.exp(sigdt*g+der)
      Se=Se*(1+sigdt*g+rdt)
      Sm=Sm*(1+sigdt*g+(sigdt*g)**2/2+der)
      
      
      
#contribution des trajectoires de S au Put
  paymc=(K-S)*(K>S)
#contribution des trajectoires du schéma d'Euler au Put
  payeul=(K-Se)*(K>Se)
  puteul=sum(payeul)
  careul=sum(payeul**2)
  ############################################################
  ## À compléter en stockant dans contputeul (resp. contcareul)
  ## la somme (resp. somme des carrés) des différences entre payoff 
  ## avec Se et payoff avec S
  ############################################################
  contputeul=sum(payeul-paymc)
  contcareul=sum((payeul-paymc)**2)
  
  #contribution des trajectoires du schéma de Milstein au Put
  #######################################################################
  ## À reprendre en remplaçant le schéma d'Euler par le schéma de Milstein
  #########################################################################
  paymil=(K-Sm)*(K>Sm)
  putmil=sum(paymil)
  carmil=sum(paymil**2)
  contputmil=sum(paymil-paymc)
  contcarmil=sum((paymil-paymc)**2)
  
  erreul[i]=np.exp(-r*T)*puteul/M-BS 
  liceul[i]=1.96*np.exp(-r*T)*np.sqrt((careul/M-(puteul/M)**2)/M)
  conterreul[i]=np.exp(-r*T)*contputeul/M
  contliceul[i]=1.96*np.exp(-r*T)*np.sqrt((contcareul/M-(contputeul/M)**2)/M)
  
  errmil[i]=np.exp(-r*T)*putmil/M-BS 
  licmil[i]=1.96*np.exp(-r*T)*np.sqrt((carmil/M-(putmil/M)**2)/M)
  conterrmil[i]=np.exp(-r*T)*contputmil/M
  contlicmil[i]=1.96*np.exp(-r*T)*np.sqrt((contcarmil/M-(contputmil/M)**2)/M)
  
  Npas[i]=N
  N=N*2#multiplication du nombre N de pas par 2
print("Erreur faible Euler")
print(erreul)
print("LIC Euler")
print(liceul)
print("Erreur faible Euler VA controle")
print(conterreul)
print("LIC Euler VA controle")
print(contliceul)

#representation graphique de (h=T/N,conterreul)
plb.plot(T/Npas,conterreul, color="r", label="erreur faible euler")
plb.plot(T/Npas,conterreul-contliceul, color="b", label="IC 95%")
plb.plot(T/Npas,conterreul+contliceul, color="b")
plb.xlabel('Pas de discretisation')
plb.legend(loc="best")
plb.show()   

print("Erreur faible Milstein")
print(errmil)
print("LIC Milstein")
print(licmil)
print("Erreur faible Milstein VA controle")
print(conterrmil)
print("LIC Milstein VA controle")
print(contlicmil)

#representation graphique de (h=T/N,conterreul)
plb.plot(T/Npas,conterrmil, color="r", label="erreur faible Milstein")
plb.plot(T/Npas,conterrmil-contlicmil, color="b", label="IC 95%")
plb.plot(T/Npas,conterrmil+contlicmil, color="b")
plb.xlabel('Pas de discretisation')
plb.legend(loc="best")
plb.show()