#c21.txt: Brownian motion done right
 Help:=proc(): print(`Bt(t,N), ExpP(F,B,t,N,M,t0)`): 
  print(`Ta(B,N,a)  , Ref(B,N,T)`):
  end:
  
  ra:=proc(): 2*rand(0..1)()-1: end:  

  #Bt(t,N): A sample Brownian path of duration t
  #with N big, dt=1/N, dx=sqrt(dt) (N should a perfect square)
  Bt:=proc(t,N) local dt,dx,i,j,P:
  dt:=1/N:
  dx:=sqrt(dt):
  P:=[seq(ra(),i=1..t*N)]:
  P:=[seq(add(P[j],j=1..i),i=1..nops(P))]:
  P:=[seq([i*dt,P[i]*dx],i=1..nops(P))]:


  end:
  
  #ExpP(F,B,t,N,M,t0): apprximates the Expected value of the
  #Stochastic process F(B,t), where F is any function
  #of B and t ,  using dt=1/N, at time t0 by simulation
  #M times
  ExpP:=proc(F,B,t,N,M,t0) local i:

   evalf(add(subs({t=t0,B=Bt(t0,N)[t0*N][2]}  ,F) , i=1..M)/M):
  end:
  
 #Ta(B,N,a): inputs a Brownian sample path, N a perfect square
 #dt=1/N, and a "real" number a, and outputs the FIRST time
 #the value of B is a 
 Ta:=proc(B,N,a) local dt,dx,i,a1:
 dt:=1/N: dx:=sqrt(dt):  
 a1:=trunc(a/dx)*dx:

  for i from 1 to nops(B) while B[i][2]<>a1 do od:
 
 if i=nops(B)+1 then
      FAIL:
   else
     i*dt:
  fi:
 
   
 end:  

  #Ref(B,N,T): inputs a sample Brownian path with dt=1/N
  #and a time, outputs the new sample path, obtained
  #by reflecting the part after T with respect to x=B[T][2]
   #WT-(B[j][2]-WT)
  Ref:=proc(B,N,T) local dx,dt,i,WT,j:
  dt:=1/N: 
  dx:=sqrt(dt):
   i:=trunc(T/dt):
   WT:=B[i][2]:
   [op(1..i,B), seq([B[j][1], 2*WT-B[j][2]], j=i+1..nops(B))]:
 end: