#Estimation of several integrals related to the Brownian motion


import numpy as np
import numpy.random as npr
import scipy.stats as sps
import matplotlib.pyplot as plt

T=2
N=int(1e4)
deltat=T/float(N)

N_LGN=int(2e4)
t=np.linspace(start=0,stop=T,num=N)


Sample=[]

for i in range(N_LGN):
    dB=np.sqrt(deltat)*npr.randn(N)    
    dX=t*dB
    XT=np.sum(dX)
    Sample.append(XT)
    
MoySample=np.mean(Sample)
StdSample=np.std(Sample)
error=1.96*StdSample*N_LGN**(-.5)

print("\n"+"X_T :") 
print("Moyenne = "+str(MoySample))
print("Moyenne theorique = "+str(0))
print("Erreur a 95% = "+str(error))

Sample2=np.array(Sample)**2
MoySample2=np.mean(Sample2)
StdSample2=np.std(Sample2)
error2=1.96*StdSample2*N_LGN**(-.5)

print("\n"+"X_T^2 :") 
print("Moyenne = "+str(MoySample2))
print("Moyenne theorique = "+str(T**3/3.))
print("Erreur a 95% = "+str(error2))