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

#Read procedures from class
read(`C16.txt`):

#Help procedure
Help:=proc():
    print(` PO(H, r, dt, C) `):
end:

##PROBLEM 1##

#Fixed version of FL (using evalf on all the exponentials)
FL:=proc(S0, Sd, Su, r, dt, C) local q, P, phi, psi, expval:
    if r <> 0 then
        expval:=evalf(exp(r*dt)):
    else
        expval:=exp(r*dt):
    fi:
    
    q:=(expval*S0-Sd)/(Su-Sd):
    
    if q < 0 or q > 1 then
        return FAIL:
    fi:
    
    P:=q*C[2]/expval + (1-q)*C[1]/expval:
    
    phi:=(C[2]-C[1])/(Su-Sd):
    psi:=P-phi*S0:
    
    return [q, P, phi, psi]:
end:

#PO(H, r, dt, C): inputs all possible stock-price histories, where
#H is a list of size, k+1, say, where H[i+1][j+1] is the price of
#the stock at time t=i if the history of downs and ups (starting
#at time 0) is coded by the binary representation of the integer
#j (j is between 0 meaning it was always down before time t=i, and
#2^i-1, meaning it was always up), r is interest rate and dt is
#the duration of one time-unit. It also inputs the list of length
#2^k-1 , where C[j+1] is the value of the claim if the history up
#to time t=k was coded by the integer j [j between 0 and 2^k-1,
#of course].
#It outputs a list FC, analogous to H where FC[k+1] is C, and
#FC[i+1][j+1] is the value of the option at time t=i if the
#history, up to that time is coded by j. In particular, FC[1]
#should be a list of length 1, whose value is the fair price of the option.
PO:=proc(H, r, dt, C) local FC, i, j, k:
    k:=nops(H)-1:
    FC[k+1]:=C:
    for i from k-1 by -1 to 0 do
        for j from 0 to 2^i-1 do
            FC[i+1][j+1]:=FL(H[i+1][j+1], H[i+2][2*j+1],
                            H[i+2][2*j+2], r, dt,
                            [FC[i+2][2*j+1], FC[i+2][2*j+2]])[2]:
        od:
    od:
    return [seq([seq(FC[i+1][j+1], j=0..2^i-1)], i=0..k)]:
end: