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

#Read packages
with(`combinat`):

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

#Help procedure
Help:=proc():
    print(` Arbi1(S0, K, r, Sd, Su, T, P, x2) `):
    print(` Goodx2(S0, K, r, Sd, Su, T, P, h) `):
    print(` EmpiOptionPrice(S0, K, r, Sd, Su, T, h1, h2) `):
end:

##PROBLEM 1##

#Arbi1(S0, K, r, Sd, Su, T, P, x2): Modification of Arbi that does
#not input x1, but inside the program x1 is set to P-x2*S0. In
#other words, the seller of the option uses the full amount he or
#she got from the buyer of the option to build his or her portfolio.
Arbi1:=proc(S0, K, r, Sd, Su, T, P, x2)
    return Arbi(S0, K, r, Sd, Su, T, P, P-x2*S0, x2):
end:

##PROBLEM 2##

#Goodx2(S0, K, r, Sd, Su, T, P, h): inputs the same as Arbi1
#(except for x2) and a mesh-size h, and by trying out all x2 from
#0 to 1, in increments of h, finds the set of x2's (that happens
#to be an interval), that allows arbitrage. Of course, it could be
#the empty set (interval).

Goodx2:=proc(S0, K, r, Sd, Su, T, P, h) local x2, S:
    S:={}:
    for x2 from 0 by h to 1 do
        if Arbi1(S0, K, r, Sd, Su, T, P, x2)[1] then
            S:=S union {x2}:
        fi:
    od:
    return S:
end:

#Goodx2(2500, 3000, 0, 2000, 4000, 1, 500, 0.001) returned the
#interval from .334 to 1.000 (presumably meaning the truth is
#(1/3, 1]
#Goodx2(2500, 3000, 0, 2000, 4000, 1, 450, 0.001) returned the
#interval from .367 to 0.900
#Goodx2(2500, 3000, 0, 2000, 4000, 1, 400, 0.001) returned the
#interval from .400 to 0.800
#Goodx2(2500, 3000, 0, 2000, 4000, 1, 200, 0.001) returned the
#empty interval

##PROBLEM 3##

#EmpiOptionPrice(S0, K, r, Sd, Su, T, h1, h2): inputs everything as
#Goodx2, except for P, and mesh sizes h1, h2, tries out all the
#P's, starting at 0, with increment h2, examines whether
#Goodx2(S0, K, r, Sd, Su, T, P, h1) is empty, and stops as soon as
#it is not empty, and returns the largest P for which it is empty.
EmpiOptionPrice:=proc(S0, K, r, Sd, Su, T, h1, h2) local P:
    for P from 0 by h2 to infinity do
        if nops(Goodx2(S0, K, r, Sd, Su, T, P, h1)) > 0 then
            return P - h2:
        fi:
    od:
end:

#EmpiOptionPrice(2500, 3000, 0, 2000, 4000, 1, 0.001, 1)
#returned 250
#EmpiOptionPrice(2500, 3000, 0.1, 2000, 4000, 1, 0.001, 1)
#returned 345
#EmpiOptionPrice(2500, 3000, 0.2, 2000, 4000, 1, 0.001, 1)
#returned 431
#EmpiOptionPrice(3500, 3000, 0, 1000, 5000, 1, 0.001, 1)
#returned 1250

##PROBLEM 4##

#See hw11.pdf