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

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

#Help procedure
Help:=proc():
    print(` CheckNewMartA(k, M, U, K) , CheckNewMart(k, M, U, K, HowMany) `):
    print(` Findphi(M, X) , BSDCP(u, d, r, dt, N, K, call) `):
    print(` BSdiscreteCall(u, d, r, dt, N, K) , BSdiscretePut(u, d, r, dt, N, K) `):
end:

##PROBLEM 1##

#CheckNewMartA(k, M, U, K): first creates a random martingale,
#using, RandMar(k, M, U, K), then generates a random phi (a random
#list of numbers with the same length as the number of non-leaves
#of that martingale), then uses procedure NewMart with a random Z0,
#and then uses, procedure IsMartingale(P) to see if the new
#martingale is indeed a martingale, as promised.
#
#NOTE: Every element of phi will be between 1 and U, inclusive
#NOTE: Z0 is between 0 and U, inclusive
CheckNewMartA:=proc(k, M, U, K) local P, ra, phi, i, Z0:
    P:=RandMar(k, M, U, K):
    ra:=rand(1..U):
    phi:=[seq(ra(), i=1..nops(P[2]))]:
    Z0:=rand(0..U)():
    P:=NewMart(P, phi, Z0):
    return IsMartingale(P):
end:

##PROBLEM 2##

#CheckNewMart(k, M, U, K, HowMany): runs CheckNewMartA(k, M, U, K)
#HowMany times, and returns true if all of them are true, and
#false otherwise.
CheckNewMart:=proc(k, M, U, K, HowMany) local i:
    return convert([seq(CheckNewMartA(k, M, U, K), i=1..HowMany)], `and`):
end:

#CheckNewMart(10,3,3,5,1000) took a really long time, but changing
#1000 to 10 made it return true.  Also, changing 10 to 3 made it
#return true.

##PROBLEM 3##

#Findphi(M, X): inputs two martingales, M, and X (with the same
#underlying trees and prob. distribution (i.e. the same sample
#space and sequence of sigma fields), and outputs Z0 and phi.
#
#Return value is Z0, phi
#
#phi may have entries of the form phi[i] for some i (instead of a
#number).  This is because these systems are often underspecified.
#Such value of phi can be anything and they won't change the
#resulting martingale.  This happens if and only if M has a node
#with exactly one child.
Findphi:=proc(M, X) local eq, var, phi, path, n, Z0, P:
    eq:={}:
    for n from 0 to M[1]-1 do
        path:=[n]:
        P:=ParentList(M):
        while path[1] <> 0 do
            path:=[P[path[1]], op(path)]:
        od:
        
        eq:=eq union {X[5][n+1] =
            Z0 + add(phi[path[j+1]+1]*(M[5][path[j+2]+1]
               - M[5][path[j+1]+1]), j=0..nops(path)-2)}:
    od:
    var:={Z0, seq(phi[j], j=1..nops(M[2]))}:
    var:=solve(eq, var):
    return op(subs(var, [Z0, [seq(phi[j], j=1..nops(M[2]))]])):
end:

FindPhi:=Findphi: #stupid capitalization conventions...

##PROBLEM 4##
#Note: the problem said "the fair price of stock" originally.  This
#should have been some sort of option, so I went with European,
#though this doesn't use any of the new martingale stuff
#This uses some stuff from hw17
#Also, no time step was given, so we added it

#Auxiliary procedure
BSDCP:=proc(u, d, r, dt, N, K, call) local L, i, j:
    L:=[seq([seq(u^i*d^(j-i), i=0..j)], j=0..N)]:
    return PogCP(L, r, dt, K, call)[1][1][2]:
end:

#BSdiscreteCall(u, d, r, N, K): computes, NUMERICALLY, the fair
#price of a European (I hope!) call option
#whose initial price is 1 dollar, lasts N
#generations, and follows the BINOMIAL model, where the price of
#the stock, from one generation to the next, may either go up by a
#factor u (u>1), or down with a factor d (d<1), with interest rate
#r, and strike price K.
BSdiscreteCall:=proc(u, d, r, dt, N, K)
    return BSDCP(u, d, r, dt, N, K, 1):
end:

#BSdiscretePut(u, d, r, N, K): computes, NUMERICALLY, the fair
#price of a European (I hope!) put option
#whose initial price is 1 dollar, lasts N
#generations, and follows the BINOMIAL model, where the price of
#the stock, from one generation to the next, may either go up by a
#factor u (u>1), or down with a factor d (d<1), with interest rate
#r, and strike price K.
BSdiscretePut:=proc(u, d, r, dt, N, K)
    return BSDCP(u, d, r, dt, N, K, -1):
end:

#BSdiscreteCall(2,1/2,0,dt,30,1.1) returned 0.9358432232
#BSdiscretePut(2,1/2,0,dt,30,1.1) returned 1.035843224