#April 16, 2015; C23.txt, Ito's famous SDE
Help:=proc(): print(`Ito(f,t,Wt), ExpRec(f,t,Wt) `):
print(`ExpRecDE(f,t,Wt,g), MakeSDE(f,t,Wt,Zt)`):
print(`Ints(f,s,B,B0,t0,N), IntW(f,s,B,B0,t0,N), Bt(t,N) `):
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:=evalf([seq([i*dt,P[i]*dx],i=1..nops(P))]):
end:

#Ito(f,t,Wt): Ito's generic SDE for an expression f of t and Wt
#[dWtComponent, dtComp]
Ito:=proc(f,t,Wt):
simplify([diff(f,Wt),diff(f,t)+1/2*diff(f,Wt,Wt)]):
end:

#expressing the expectation of f(t,Wt) in terms of
#hopefully simpler things
ExpRec:=proc(f,t,Wt) local sid:
option remember:

if diff(f,Wt)=0 then 
 RETURN(f):
fi:

if diff(f/Wt,Wt)=0 then
 RETURN(0):
fi:

sid:=Ito(f,t,Wt)[2]:
int(ExpRec(sid, t,Wt),t):

end:

#MakeSDE(f,t,Wt): inputs an expression f in t and Wt
#and outputs an SDE (Eq. 4.14) outputs the
#pair [sig(t,Zt),mu(t,Zt)]
MakeSDE:=proc(f,t,Wt,Zt) local ito,Wt0:

ito:=Ito(f,t,Wt):

Wt0:=solve(f=Zt,Wt):

simplify(subs(Wt=Wt0,ito)):


end:


#written by Anthony Zaleski
#Ints(f,s,B,B0,t0,N): inputs an expression f in s and B; B0, a specific discrete 
#sample Brownian motion (gotten from Bt(t0,N); a time t0; and a parameter N 
#(a perfect square) and outputs a discrete approximation to the usual integral 
#of f(s,B_0) with respect to s from 0 to t0.

Ints:=proc(f,s,B,B0,t0,N) local i:
if t0>B0[-1][1] then return FAIL: fi:
evalf(add(subs({s=B0[i][1],B=B0[i][2]},f)/N , i=1..N*t0)):
end:


#IntW(f,s,B,B0,t0,N): inputs an expression, f in in s and B; B0, a specific 
#discrete sample Brownian motion (gotten from Bt(t0,N)); a time t0; and a 
#parameter N (a perfect square) and outputs a discrete approximation to the 
#STOCHASTIC integral of f(s,B_0) with respect to dW_s, from 0 to t0.

IntW:=proc(f,s,B,B0,t0,N) local i:
if t0>B0[-1][1] then return FAIL: fi:
evalf(add(subs({s=B0[i][1],B=B0[i][2]},f)*(B0[i][2]-B0[i-1][2]) , i=2..N*t0)):
end:
#End written by Anthony Zaleski