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

#Read packages
with(`LinearAlgebra`):

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

#Help procedure
Help:=proc():
    print(` PrC(d, u, k) , PrP(d, u, k) `):
    print(` ECB(d, u, k, q) , ECS(d, u, k, q) , EPB(d, u, k, q) , EPS(d, u, k, q) `):
    print(` Problem3(part, d, u, k, q, e) `):
    print(` PriceEO(CER, FERD, FERU, u, r, T, Pou, Dol) `):
    print(` LogNormalStdMoments(mu, sigma, k) `):
end:

##PROBLEM 1##

#For the sake of simplicity, we assume that r=0, and that the
#current stock price is 1.

#PrC(d, u, k): the price of a European Call option with strike price k
PrC:=proc(d, u, k)
    return (u-k)*(1-d)/(u-d):
end:

#PrP(d, u, k): the price of a European Put option with strike price k
PrP:=proc(d, u, k)
    return (k-d)*(u-1)/(u-d):
end:

#PROBLEM 2##

#The buyer (resp. seller) believes that the probability that the
#stock will go up is q

#Expected Buyer Call Gain
ECB:=proc(d, u, k, q)
    return q*(u-k)-PrC(d, u, k):
end:

#Expected Seller Call Gain
ECS:=proc(d, u, k, q)
    return PrC(d, u, k)-q*(u-k):
end:

#Expected Buyer Put Gain
EPB:=proc(d, u, k, q)
    return (1-q)*(k-d)-PrP(d, u, k):
end:

#Expected Seller Put Gain
EPS:=proc(d, u, k, q)
    return PrP(d, u, k)-(1-q)*(k-d):
end:

##PROBLEM 3##

Problem3:=proc(part, d, u, k, q, e)
    if part = `a` then
        return solve(ECB(d, u, k, q) + ECS(d, u, k+e, q) > 0, q) assuming d<k, u>k+e, e>0, d>0:
    elif part = `b` then
        return solve(ECB(d, u, k+e, q) + ECS(d, u, k, q) > 0, q) assuming d<k, u>k+e, e>0, d>0:
    elif part = `c` then
        return solve(ECB(d, u, k, q) + 2*EPB(d, u, k, q) > 0, q) assuming d<k, u>k+e, d>0:
    elif part = `d` then
        return solve(2*ECB(d, u, k, q) + EPB(d, u, k, q) > 0, q) assuming d<k, u>k+e, d>0:
    elif part = `e` then
        return [(solve(ECB(d, u, k, q) + EPB(d, u, k+e, q) > 0, q) assuming d<k, u>k+e, d>0), (solve(ECB(d, u, k+e, q) + EPB(d, u, k, q) > 0, q) assuming d<k, u>k+e, d>0)]:
    else
        return FAIL:
    fi:
end:

#In all cases, if there are two strike prices, k is the lower one
#a: Gains when q>(d-1)/(d-u)
#b: Gains when q<(d-1)/(d-u)
#c: If k>(u+2*d)/3, then gains when q<(d-1)/(d-u),
### otherwise gains when q>(d-1)/(d-u)
#d: If k>(u+2*d)/3, then gains when q>(d-1)/(d-u),
### otherwise gains when q<(d-1)/(d-u)
#e: If the put strike is higher and k>(u+d-e)/2, then gains when q<(d-1)/(d-u)
### otherwise if the put strike is higher, then gains when q>(d-1)/(d-u)
### If the call strike is higher and k>(u+d-e)/2, then gains when q<(d-1)/(d-u)
### otherwise if the call strike is higher, then gains when q>(d-1)/(d-u)

##PROBLEM 4##

#PriceEO(CER, FERD, FERU, u, r, T, Pou, Dol): inputs
#CER: the current Exchange Rate between a Dollar and a Pound,
#FERD: the future exchange rate in case it is going to go down
#FERU: the future exchange rate in case it is going to go up
#u: the interest rate in the UK
#r: the interest rate in the USA
#T: the time for exercising
#Pou: an amount in Pounds
#Dol: an amount in Dollars
#
#and outputs the fair price of an option to buy, at time T, Pou
#Pounds for Dol Dollars.
PriceEO:=proc(CER, FERD, FERU, u, r, T, Pou, Dol) local p:
    p:=solve(CER*exp(-u*T)=exp(-r*T)*(FERU*p+FERD*(1-p)), p):
    return exp(-r*T)*p*(FERU*Pou-Dol):
end:

##PROBLEM 5##

lognormal:=(mu, sigma, x)->(1/(x*sigma*sqrt(2*Pi)))*exp(-(log(x)-mu)^2/(2*sigma^2)):

#The answer to Problem 5 is the first moment returned here
#return the first k moments
LogNormalStdMoments:=proc(mu, sigma, k) local L, i, x:
    L:=[int(x*lognormal(mu, sigma, x), x=0..infinity) assuming sigma>0]:
    L:=[op(L), (int(x^2*lognormal(mu, sigma, x), x=0..infinity) assuming sigma>0) - L[1]^2]:
    for i from 3 to k do
        L:=[op(L), int(lognormal(mu, sigma, x)*(x-L[1])^i, x=0..infinity) assuming sigma>0]:
    od:
    return L:
end:

#The moments are
#[exp(1/2*sigma^2+mu),
#exp(2*mu+2*sigma^2)-exp(sigma^2+2*mu),
#exp(3*mu+9/2*sigma^2)-3*exp(5/2*sigma^2+3*mu)+2*exp(3*mu+3/2*sigma^2),
#exp(4*mu+8*sigma^2)-4*exp(5*sigma^2+4*mu)+6*exp(3*sigma^2+4*mu)-3*exp(4*mu+2*sigma^2)]