import numpy as np
import scipy.stats as st

np.random.seed(42)

m = 1000
r = 0.5
sigma = 0.3
T = 2
S0 = 5
K = 2*S0
Blist = np.arange(0,S0)
nB = np.size(Blist)
dt = T/m
t=np.arange(0,m)*dt
N = np.int64(1E5)
M = np.int64(N/2)
#Defines for the confidence intervals
alpha=0.05

cval = st.norm.ppf(1-alpha/2);

#Defines the gen psi fucntion 
def genPsi(xi,t,r,sigma,K,B,S0):
    dt=np.diff(t)
    W=np.append([0],np.cumsum(np.sqrt(dt)*xi))
    S=S0*np.exp((r-sigma**2/2)*t+sigma*W);
    rval = (abs(S[-1]-K)+ (S[-1]-K))/2*(B<=min(S));
    return rval

#Begins the for loop
for j in range(nB):
    B=Blist[j]
    print('B =  ',B)
    psi = np.zeros(N)
    for i in range(N):
        xi = np.random.standard_normal(m-1)
        psi[i] = genPsi(xi,t,r,sigma,K,B,S0)
    print("MC confidence interval (alpha= ",alpha, ' N= ',N,') mean= ',
          np.mean(psi), " +- ", cval*np.std(psi)/np.sqrt(N))
    NvarMC = np.var(psi)
    print("N*variance MC estimator: " ,NvarMC)
    psiAV = np.zeros(M)
    for i in range(M):
        xi = np.random.standard_normal(m-1);       
        psiAV[i] = (genPsi(xi,t,r,sigma,K,B,S0) + genPsi(-xi,t,r,sigma,K,B,S0))/2
    print("AV confidence interval (alpha= ",alpha, ' N/2= ',M,') mean= ',
          np.mean(psiAV), " +- ", cval*np.std(psiAV)/np.sqrt(M))
    NvarAV = 2*np.var(psiAV)
    print('N*variance AV estimator: ',NvarAV)