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

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

#Help procedure
Help:=proc():
    print(` Tab(B, N, a, b) , Check335(a, b, t, N, theta, M) `):
    print(` Check334(a, x, t0, N, M) `):
end:

##PROBLEM 1##

#Tab(B, N, a, b): finds the first time B(t, N) meets the line
#x=a+b*t.
#
#Note: "Meets" means "Crosses"
#It does not need to intersect exactly
#
#Note: may not return the same thing as Ta(B, N, a, 0), due to
#truncation in the former procedure
Tab:=proc(B, N, a, b) local dt, i, sgn:
    if a = 0 then
        return 0:
    fi:
    sgn:=sign(a):
    dt:=1/N:
    #dx:=sqrt(dt):
    
    for i from 1 to nops(B) do
        if sgn*B[i][2] >= sgn*(a+b*dt*i) then
            return i*dt:
        fi:
    od:
    return FAIL:
end:

##PROBLEM 2##

#Check335(a, b, t, N, theta, M): verifies Prop. 3.3.5 on p. 61 of
#Alison Etheridge's book (approximately) by picking
#M random Bt(t, N)'s.
#Return value is the calculated expectation, the
#expected expectation, and their difference
#
#NOTE: I added an argument t to this procedure
Check335:=proc(a, b, t, N, theta, M) local i, T, B, val, exval:
    val:=0:
    for i from 1 to M do
        B:=Bt(t, N):
        T:=Tab(B, N, a, b):
        if T <> FAIL then
            val:=val + exp(-theta*T):
        fi:
    od:
    val:=evalf(val/M):
    exval:=evalf(exp(-a*(b+sqrt(b^2+2*theta)))):
    return val, exval, abs(val-exval):
end:

##PROBLEM 3##

#Check334(a, x, t0, N, M): verifies Prop. 3.3.4 on p. 60 of
#Alison Etheridge's book (approximately) by picking
#M random Bt(t, N)'s.
#Return value is the calculated probability, the
#expected probability, and their difference
Check334:=proc(a, x, t0, N, M) local i, val, t, Mt, B, exval:
    val:=0:
    for i from 1 to M do
        B:=Bt(t0, N):
        Mt:=max(seq(B[t][2], t=1..nops(B))):
        if Mt >= a and B[-1][2] <= x then
            val:=val+1:
        fi:
    od:
    val:=evalf(val/M):
    exval:=evalf(1-PHI((2*a-x)/sqrt(t0))):
    return val, exval, abs(val-exval):
end: