# Ej 19

import random
import numpy as np
import plotly.express as px
import matplotlib.pyplot as plt
import scipy as sc
import scipy.stats as stats


def simulate_1d_bm(nsteps=1000, t=0.001):
    steps = [ np.random.randn()*np.sqrt(t) for i in range(nsteps) ]
    steps[0] = 0
    y = np.cumsum(steps)
    x = [ t*i for i in range(nsteps) ]
    return x, y

def geometric_bm(nsteps=1000, t=0.001,S_0=0, sigma=1,r=0):
    steps = [ np.random.randn()*np.sqrt(t) for i in range(nsteps) ]
    steps[0] = 0
    aux = np.cumsum(steps)
    y= np.zeros_like(aux)
    mu = r - sigma**2/2
    for k in range(len(aux)):
        y[k] = S_0*np.exp(sigma*aux[k] + k*t*mu)
    x = [ t*i for i in range(nsteps) ]
    return x, y


def integrate_gm_b(F,dt,t_final):
    return np.sum(F)*dt/t_final
#%% Example

S_0=1
sigma=1 #1 
r=0.08
fig = plt.figure()
fig.set_figheight(8)
fig.set_figwidth(12)
nsims = 5

dt = 0.001
n_steps = 2000

for i in range(nsims):
    x, y = geometric_bm(n_steps,dt,S_0,sigma,r)
    z, w = simulate_1d_bm(n_steps,dt)
    print(integrate_gm_b(x,dt,2))
    plt.plot(x,y,'k'),plt.plot(z,w+1,'b')

plt.show()

#%%

n_sims=100
S_0=1
K= 1
sigma=0.4 #1 
r=0.08
dt = 0.001
T = 1
n_steps = 1000

fig = plt.figure()
fig.set_figheight(8)
fig.set_figwidth(12)

fun = np.zeros(n_sims)

for i in range(n_sims):
    x, y = geometric_bm(n_steps,dt,S_0,sigma,r)
    integral = integrate_gm_b(x,dt,T)
    
    if (integral -K)> 0:
        fun[i] = integral -K
    
    
E_St_K = np.mean(fun)