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

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

#Help procedure
Help:=proc():
    print(` MyTaylor(f, x, a, N) , Taylor2Start(f, x, y, h, k, N) `):
    print(` ExpWtr(t, r) , Ints(f, s, B, B0, t0, N) `):
    print(` IntW(f, s, B, B0, t0, N) `):
end:

##PROBLEM 1##

#MyTaylor(f, x, a, N): inputs an expression f in the variable x, a
#number (or symbol) a, and a positive integer N, and outputs the
#polynomial of degree N consisting of the first N terms of the
#Taylor expansion of f at x=a.
MyTaylor:=proc(f, x, a, N) local n:
    return simplify(subs(x=a, f)) + add(simplify(subs(x=a, diff(f, x$n))/n!)*(x-a)^n, n=1..N-1):
end:

##PROBLEM 2##

#Taylor2Start(f, x, y, h, k, N): inputs an expression f in the
#variables x, and y, symbols h and k, and a positive integer N,
#and outputs all the terms of total degree <=N in the Taylor
#expansion of f(x+h,y+k) around the point (x,y).
Taylor2Start:=proc(f, x, y, h, k, N) local d, i:
    return f + add(add(diff(f, x$i, y$(d-i))*h^i*k^(d-i)
        *binomial(d, i)/d!, i=0..d), d=1..N):
end:

#Taylor2Start(exp(x+y+x*y), x, y, h, k, 3); outputs
#exp(x*y+x+y)+(x+1)*exp(x*y+x+y)*k+(y+1)*exp(x*y+x+y)*h+1/2*(x+1)^2*exp(x*y+x+y)*k^2+(exp(x*y+x+y)+(y+1)*(x+1)*exp(x*y+x+y))*h*k+1/2*(y+1)^2*exp(x*y+x+y)*h^2+1/6*(x+1)^3*exp(x*y+x+y)*k^3+1/2*(2*(x+1)*exp(x*y+x+y)+(y+1)*(x+1)^2*exp(x*y+x+y))*h*k^2+1/2*(2*(y+1)*exp(x*y+x+y)+(y+1)^2*(x+1)*exp(x*y+x+y))*h^2*k+1/6*(y+1)^3*exp(x*y+x+y)*h^3
#Taylor2Start(f(x,y), x, y, h, k, 3); outputs
#f(x,y)+diff(f(x,y),y)*k+diff(f(x,y),x)*h+1/2*diff(diff(f(x,y),y),y)*k^2+diff(diff(f(x,y),x),y)*h*k+1/2*diff(diff(f(x,y),x),x)*h^2+1/6*diff(diff(diff(f(x,y),y),y),y)*k^3+1/2*diff(diff(diff(f(x,y),x),y),y)*h*k^2+1/2*diff(diff(diff(f(x,y),x),x),y)*h^2*k+1/6*diff(diff(diff(f(x,y),x),x),x)*h^3

##PROBLEM 3##

#ExpWtr(t, r): uses Ito's formula, and a recursion (generalizing
#example 4.3.2 on p. 87) to compute the Expectation of the
#Stochastic process
#(W_t)^r
#for any specific (numeric) NON-NEGATIVE integer r, yielding an
#expression in t.
ExpWtr:=proc(t, r) local s:
    if r = 0 then
        return 1:
    elif r = 1 then
        return 0:
    else
        return 1/2*r*(r-1)*int(ExpWtr(s, r-2), s=0..t):
    fi:
end:

#Conjecture: ExpWtr(t, r)=0 if r is odd and =r!!*t^(r/2) if r is even
#(where r!!=product of odd positive integers less than r)
#Proof: The zero case is easy by induction.  The other case follows
#as well by induction.  If r=0, W_t^0=1 always. Later, we have
#ExpWtr(t, r)=1/2*r*(r-1)*int((r-2)!!*s^(r/2-1), s=0..t)
#=1/2*r*(r-1)*(r-2)!!*t^(r/2)*(1/(r/2))
#=1/2*(r-1)*(r-2)!!*t^(r/2)
#=1/2*r!!*t^(r/2), as required.

##PROBLEM 4##

#Ints(f, s, B, B0, t0, N): inputs an expression, f in s and B, B0,
#a specific discete 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 (classical) integral of
#f(s, B_0) with respect to s from 0 to t0.
Ints:=proc(f, s, B, B0, t0, N) local n, dt:
    dt:=1/N:
    return subs({s=0, B=0}, f) + add(subs({s=n*dt, B=B0[n][2]}, f)*dt, n=1..N*t0):
end:

##PROBLEM 5##

#IntW(f, s, B, B0, t0, N): inputs an expression, f in s and B, B0,
#a specific discete 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 n, dt:
    dt:=1/N:
    return subs({s=dt, B=B0[1][2]}, f)*B0[1][2] +
        add(subs({s=B0[n][1], B=B0[n][2]}, f)*(B0[n][2]-B0[n-1][2]), n=2..N*t0):
end: