#Estimation of several integrals or distributions 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=5.

N=int(2e4)
deltat=T/float(N)
N_LGN=int(2e4)
t=np.linspace(start=0,stop=T,num=N)
expmt=np.exp(-t)

Sample=[]

for i in range(N_LGN):
    dB=np.sqrt(deltat)*npr.randn(N)
    Sample.append(np.sum(expmt*dB))

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)

sigma2=(1-np.exp(-2*T))/2.

print("\n"+"X_T^2 :") 
print("Moyenne = "+str(MoySample2))
print("Moyenne theorique = "+str(sigma2))
print("Erreur a 95% = "+str(error2))

Sample=sigma2**(-.5)*np.array(Sample)
plt.figure()
plt.hist(Sample, bins=round(4*len(Sample)**.3),normed=1,color='b',histtype='step')
x=np.linspace(-4,4,1000)
y=sps.norm.pdf(x)
plt.plot(x,y,color='r')