#Nathan Fox
#Homework 23
#I give permission for this work to be posted online

#Read procedures from class/last homework
read(`C23.txt`):
read(`hw22.txt`):

#Help procedure
Help:=proc():
    print(` ExpftW(f, t, Wt, t0, N, M) , ExpInts(f, s, B, t0, N, M) `):
    print(` ExpIntW(f, s, B, t0, N, M) , CheckIto(f, s, Wt, t0, N, M) `):
    print(` ExpftWexact(f, t, Wt) `):
end:

##PROBLEM 1##

#ExpftW(f, t, Wt, t0, N, M): inputs an expression f in t and Wt,
#a real number t0, a large integer N (a perfect square), and a
#fairly large integer M, and empirically estimates the expectation
#of the value of the Stochastic Process f(t, W_t) by taking M
#random discrete Brownian motions (using Bt(t0, N) M times),
#and averaging it.
ExpftW:=proc(f, t, Wt, t0, N, M) local i:
    return add(subs({t=t0, Wt=Bt(t0, N)[-1][2]}, f), i=1..M)/M:
end:

##PROBLEM 2##

#ExpInts(f, s, B, t0, N, M): estimates the classical integral
#Int(f(s,B)ds) using Ints(f, s, B, B0, t0, N) with M
#randomly-picked B0's from Bt(t0, N)
#
#NOTE: There is no reason for B0 to be an argument to this procedure
#so I made it not be one
ExpInts:=proc(f, s, B, t0, N, M) local i:
    return add(Ints(f, s, B, Bt(t0, N), t0, N), i=1..M)/M:
end:

##PROBLEM 3##

#ExpIntW(f, s, B, t0, N, M): estimates the stochastic integral
#Int(f(s,W)dW) using IntW(f, s, B, B0, t0, N) with M
#randomly-picked B0's from Bt(t0, N)
#
#NOTE: There is no reason for B0 to be an argument to this procedure
#so I made it not be one
ExpIntW:=proc(f, s, B, t0, N, M)
    return add(IntW(f, s, B, Bt(t0, N), t0, N), i=1..M)/M:
end:

##PROBLEM 4##

#CheckIto(f, s, Wt, B0, t0, N, M): empirically verifies Ito's SDE
#for the Stochastic process f(t, W)
#
#The two numbers returned are the two sides of Ito's SDE
#They should be close to each other
#
#NOTE: There is no reason for B0 to be an argument to this procedure
#so I made it not be one
CheckIto:=proc(f, s, Wt, t0, N, M)
    return [ExpftW(f, s, Wt, t0, N, M) - subs({s=0, Wt=0}, f),
        ExpIntW(diff(f, Wt), s, Wt, t0, N, M) + ExpInts(diff(f, s), s, Wt, t0, N, M) + 1/2*ExpInts(diff(f, Wt$2), s, Wt, t0, N, M)]:
end:

#This does not appear to be working, but it might just be because
#everything converges so slowly

##PROBLEM 5##

#simplify(Ito(f(t, S0*exp(nu*t+sigma*Wt)), t, Wt)); returned
#[D[2](f)(t,S0*exp(Wt*sigma+nu*t))*S0*sigma*exp(Wt*sigma+nu*t), D[1](f)(t,S0*exp(Wt*sigma+nu*t))+D[2](f)(t,S0*exp(Wt*sigma+nu*t))*S0*nu*exp(Wt*sigma+nu*t)+1/2*D[2,2](f)(t,S0*exp(Wt*sigma+nu*t))*S0^2*sigma^2*exp(2*Wt*sigma+2*nu*t)+1/2*D[2](f)(t,S0*exp(Wt*sigma+nu*t))*S0*sigma^2*exp(Wt*sigma+nu*t)],
#which looks a lot like what we were trying to prove on p. 89
#if we replace each instance of S0*exp(Wt*sigma+nu*t) with St, yielding
#[D[2](f)(t,St)*St, D[1](f)(t,St)+D[2](f)(t,St)*St+1/2*D[2,2](f)(t,St)*sigma^2*St^2+1/2*D[2](f)(t,St)*sigma^2*St]
#Though, this appears to have an extra derivative everywhere.

##PROBLEM 7##

#ExpftWexact(f, t, Wt): the expected value of the Stochastic process
#f(t,Wt) is given by the following integral formula:
#Int(f(t,x)*exp(-x^2/(2*t))/sqrt(2*Pi*t),x=-infinity..infinity)
ExpftWexact:=proc(f, t, Wt)
    return int(f*exp(-Wt^2/(2*t))/sqrt(2*Pi*t), Wt=-infinity..infinity) assuming t > 0:
end:

#ExpftWexact(exp(a*Wt), t, Wt); returns exp(a*t^2/2)
#ExpftWexact(exp(nu*t+sig*Wt), t, Wt); returns exp(nu*t+1/2*t*sig^2)
#ExpftWexact(cos(a*Wt), t, Wt); returns 1/(Pi*t)^(1/2)*Pi^(1/2)*t^(1/2)*exp(-1/2*t*a^2)
#ExpftWexact(cos(t+a*Wt), t, Wt); returns 1/2*exp(-I*t-1/2*t*a^2)*(exp(2*I*t)+1)