# -*- coding: utf-8 -*-
"""
Created on Fri Apr 22 14:16:24 2022

@author: turinici
"""

import numpy as np
import matplotlib.pyplot as plt

T=1
N=500
h=T/N
M=200
trange=np.linspace(0,T,N)

W = np.zeros((N+1,M))
dW= np.random.randn(N,M)* np.sqrt(h)
W[0,:]=0
W[1:,:]=W[0,:]+ np.cumsum(dW,axis=0)
plt.figure('Brownian')
plt.plot(W)
plt.show()

S0=100
mu=0.1
sigma=0.25
K=110
r=0.05

Y=np.zeros_like(W)
Y[0,:]=S0
for ii in range(N):
    Y[ii+1,:]=Y[ii,:]+ mu*Y[ii,:]*h+ sigma*dW[ii,:]*Y[ii,:]
    
plt.figure("B&S historique")
plt.plot(Y)
plt.show()

St=np.zeros_like(W)
St[0,:]=S0
for ii in range(N):
    #Euler Maruyama
    St[ii+1,:]=St[ii,:]*(1+ r*h+ sigma*dW[ii,:])
    #Milshtein
    St[ii+1,:]=St[ii,:]*(1+ r*h+ sigma*dW[ii,:] + sigma*sigma/2*(dW[ii,:]**2-h) )
    
plt.figure("B&S historique")
plt.plot(St)
plt.show()

prix_option=np.exp(-r*T)*np.mean(np.maximum(0,St[-1,:]-K)  )
print(prix_option)