Floating call
import numpy as np
#input
s0 = 50
r = 0.06
sigma = 0.3
T = 1
M = 50
dt = T/M
u = np.exp(sigma*np.sqrt(dt))
d = 1/u
a = np.exp(r*dt)
p = (a-d)/(u-d)
#option price at time t = M
s = np.zeros(M+1)
min1 = np.zeros(M+1)
min2 = np.zeros(M-1)
for n in range(M+1):
s[n] = s0*(np.power(u,M-n)*np.power(d,n))
min1[n] = s0*np.power(d,n) # minimum from upward path
for i in range(M-1):
min2[i] = s0*np.power(d,i) #minimum from down path
s1 = np.zeros(M-1)
for i in range(M-1):
s1[i] = s[i+1] #option price at t = M excluding first and last
#value of the option
vm1 = np.zeros(M+1)
vm2 = np.zeros(M-1)
for i in range(M+1):
vm1[i] = max(s[i]-min1[i],0)
for i in range(M-1):
vm2[i] = max(s1[i]-min2[i],0)
#value at time t = M-1
x = np.zeros(M)
for i in range(M):
if i == 0:
x[i] = np.exp(-r*dt)*(p*vm1[0]+(1-p)*vm2[0])
elif i == M-1:
x[i] = np.exp(-r*dt)*(p*vm1[M-1]+(1-p)*vm1[M])
else:
x[i] = np.exp(-r*dt)*(p*vm1[i]+(1-p)*vm2[i])
#value at time t = M-1 excluding first and last
x1 = np.zeros([M-1,M])
x = np.vstack([x,x1])
for i in range(M-1):
for j in range(M-1-i):
x[i+1,j] = np.exp(-r*dt)*(p*x[i,j]+(1-p)*x[i,j+1])
print(x[M-1,0]) # result
Floating put
import numpy as np
s0 = 9
r = 0.06
sigma = 0.3
T = 3
M = 1000
dt = T/M
u = np.exp(sigma*np.sqrt(dt))
d = 1/u
a = np.exp(r*dt)
p = (a-d)/(u-d)
s = np.zeros(M+1)
max1 = np.zeros(M+1)
max2 = np.zeros(M-1)
for n in range(M+1):
s[n] = s0*(np.power(u,M-n)*np.power(d,n))
max1[n] = s0*np.power(u,M-n)
for i in range(M-1):
max2[i] = s0*np.power(u,i+(M-1)-(2*i+1))
s1 = np.zeros(M-1)
for i in range(M-1):
s1[i] = s[i+1]
#payoff
vmax1 = np.zeros(M+1)
vmax2 = np.zeros(M-1)
for i in range(M+1):
vmax1[i] = max(max1[i]-s[i],0)
for i in range(M-1):
vmax2[i] = max(max2[i]-s1[i],0)
#value at time = M-1
x = np.zeros(M)
for i in range(M):
if i == 0:
x[i] = np.exp(-r*dt)*(p*vmax1[0]+(1-p)*vmax1[1])
elif i == M-1:
x[i] = np.exp(-r*dt)*(p*vmax2[M-2]+(1-p)*vmax1[M])
else:
x[i] = np.exp(-r*dt)*(p*vmax2[i-1]+(1-p)*vmax1[i+1])
#value untill t = 0
x1 = np.zeros([M-1,M])
x = np.vstack([x,x1])
for i in range(M-1):
for j in range(M-1-i):
x[i+1,j] = np.exp(-r*dt)*(p*x[i,j]+(1-p)*x[i,j+1])
#result
print(x[M-1,0])
Fixed call
import numpy as np
# input
s0 = 9
r = 0.06
k = 10
sigma = 0.3
T = 3
M = 7
dt = T/M
u = np.exp(sigma*np.sqrt(dt))
d = 1/u
a = np.exp(r*dt)
p = (a-d)/(u-d)
# option price at time t (S(t))
s = np.zeros(M+1)
max1 = np.zeros(M+1)
max2 = np.zeros(M-1)
for n in range(M+1):
s[n] = s0*(np.power(u,M-n)*np.power(d,n))
max1[n] = s0*np.power(u,M-n)
for i in range(M-1):
max2[i] = s0*np.power(u,i+(M-1)-(2*i+1))
s1 = np.zeros(M-1)
for i in range(M-1):
s1[i] = s[i+1]
#cari max
vmax1 = np.zeros(M+1)
vmax2 = np.zeros(M-1)
for i in range(M+1):
vmax1[i] = max(max1[i]-k,0)
for i in range(M-1):
vmax2[i] = max(max2[i]-k,0)
#value at t = M-1 until t = 0
x = np.zeros(M)
for i in range(M):
if i == 0:
x[i] = np.exp(-r*dt)*(p*vmax1[0]+(1-p)*vmax1[1])
elif i == M-1:
x[i] = np.exp(-r*dt)*(p*vmax2[M-2]+(1-p)*vmax1[M])
else:
x[i] = np.exp(-r*dt)*(p*vmax2[i-1]+(1-p)*vmax1[i+1])
x1 = np.zeros([M-1,M])
x = np.vstack([x,x1])
for i in range(M-1):
for j in range(M-1-i):
x[i+1,j] = np.exp(-r*dt)*(p*x[i,j]+(1-p)*x[i,j+1])
print(x[M-1,0])
Fixed put
import numpy as np
s0 = 9
r = 0.06
k = 10
sigma = 0.3
T = 3
M = 7
dt = T/M
u = np.exp(sigma*np.sqrt(dt))
d = 1/u
a = np.exp(r*dt)
p = (a-d)/(u-d)
s = np.zeros(M+1)
min1 = np.zeros(M+1)
min2 = np.zeros(M-1)
for n in range(M+1):
s[n] = s0*(np.power(u,M-n)*np.power(d,n))
min1[n] = s0*np.power(d,n)
for i in range(M-1):
min2[i] = s0*np.power(d,i)
s1 = np.zeros(M-1)
for i in range(M-1):
s1[i] = s[i+1]
vm1 = np.zeros(M+1)
vm2 = np.zeros(M-1)
for i in range(M+1):
vm1[i] = max(k-min1[i],0)
for i in range(M-1):
vm2[i] = max(k-min2[i],0)
x = np.zeros(M)
for i in range(M):
if i == 0:
x[i] = np.exp(-r*dt)*(p*vm1[0]+(1-p)*vm2[0])
elif i == M-1:
x[i] = np.exp(-r*dt)*(p*vm1[M-1]+(1-p)*vm1[M])
else:
x[i] = np.exp(-r*dt)*(p*vm1[i]+(1-p)*vm2[i])
x1 = np.zeros([M-1,M])
x = np.vstack([x,x1])
for i in range(M-1):
for j in range(M-1-i):
x[i+1,j] = np.exp(-r*dt)*(p*x[i,j]+(1-p)*x[i,j+1])
print(x[M-1,0])