2 minute read

[Notice] Journey to the academic researcher This is the story of how I became the insightful researcher.

Implementation of the BDT’s option pricing model in python. I treated binomial tree as a network with nodes (i,j) with i representing the time steps and j representing the number of ordered price outcome (lowest - or bottom of tree - to highest).

import numpy as np
import pandas as pd
import os    # import and read files
import os.path
from scipy.optimize import minimize

# Initialise parameters
spot = [0.05, 0.0551, 0.0604] # given spot rates
N = len(spot)-1 # the number of steps
prin = 10000 # principle payment at the mature date
sigma = 0.1
u = np.exp(sigma)
d = 1/u
print(u, d)

1.1051709180756477 0.9048374180359595

def temp_forward(sp, N):
    """return a tree of forward rates with the given spot rates"""
    temp_f = np.zeros( (N+1, N+1) )
    temp_f[0,0] = sp[0]
    for i in range(1, N+1 ):
        M = i + 1
        temp_f[0, i] = u * temp_f[0, i-1]
        for j in range(1, M ):
            temp_f[j, i] = d * temp_f[j - 1, i - 1]
    
    return temp_f

temp_f = temp_forward(spot, N)
print(temp_f)

[[0.05       0.05525855 0.06107014]
 [0.         0.04524187 0.05      ]
 [0.         0.         0.04093654]]

def discount_zero(sp):
    """return the prices of bonds with the given spot rates"""
    dis_zero = []
    for idx, i in enumerate(sp):
        temp = prin*(1/(1+i)**(idx+1))
        dis_zero.append(temp)
    
    return dis_zero

discount_zero(spot)

[9523.809523809523, 8982.821171945383, 8386.694875533425]

def forward_discount(fd_rates, N):
    """calculating present value of bonds with the given forward rates and N"""
    p_forward = np.zeros( (N+1, N+1) )
    for i in range(0, N+1):
        for j in range(0, N+1):
            if fd_rates[j, i] == 0:
                continue
            else:
                p_forward[j, i] = prin/(1+fd_rates[j, i])
    return p_forward

p_forward = forward_discount(temp_f, N)
print(forward_discount(temp_f, N))

[[9523.80952381 9476.35064299 9424.44768045]
 [   0.         9567.16361867 9523.80952381]
 [   0.            0.         9606.73358871]]

def rate_gap(adj, fd_rates, N):
    """difference between the price of bonds with spot rates and the price of bonds with forward rates in order to use a solver"""
    ab_forward = np.zeros( (N+1, N+1) )
    adj_rates = np.zeros( (N+1, N+1) )
    # only last comluns should be adjusted
    for i in range(0, N+1):
        for j in range(0, N+1):
            if i == N:
                adj_rates[j, i] = fd_rates[j, i] + adj
            else:
                adj_rates[j, i] = fd_rates[j, i]

    # calculate a new prices of bonds based on adjusted forward rates
    p_forward = forward_discount(adj_rates, N)
    for i in range(0, N+1):
        ab_forward[i, N] = p_forward[i, N]

    # average each nods' value and discount
    for i in range(N-1, -1, -1):
        for j in range(N-1, -1, -1):
            if j > i:
                continue
            else:
                ab_forward[j, i] = (ab_forward[j, i+1] + ab_forward[j+1, i+1])*0.5/(1+adj_rates[j, i])
                
    target = ab_forward[0,0]
    base = discount_zero(spot)[N]
#     print("price: \n", ab_forward)
#     print("adj_rates: \n", adj_rates)
#     print("gap: ", abs(base - target))
    return abs(base - target)

adj = 0.02
N = len(spot)-1
temp_f = temp_forward(spot, N)
print(temp_f)
print(discount_zero(spot)[N])
print(rate_gap(adj, temp_f, N))

[[0.05       0.05525855 0.06107014]
 [0.         0.04524187 0.05      ]
 [0.         0.         0.04093654]]
8386.694875533425
84.98338013316788

adj = np.array([adj])
# bnds = Bounds(-0.05,0.05)
temp = minimize(rate_gap, adj, args=(temp_f, N), method="L-BFGS-B")
print(temp)
print(temp.x[0])

      fun: 2.2222817278816365e-05
 hess_inv: <1x1 LbfgsInvHessProduct with dtype=float64>
      jac: array([3312.06319915])
  message: 'ABNORMAL_TERMINATION_IN_LNSRCH'
     nfev: 116
      nit: 3
     njev: 58
   status: 2
  success: False
        x: array([0.03084632])
0.030846318327330674

temp_f = temp_forward(spot, N)
final_rates = np.zeros( (N+1, N+1) )

for i in range(0, N+1):
    for j in range(0, N+1):
        final_rates[j, i] = temp_f[j, i]

for i in range(0, N+1):
    if i == 0:
        continue
    else:
        temp_res = minimize(rate_gap, adj, args=(final_rates, i), method="L-BFGS-B")
        print(temp_res.x[0])
        for j in range(0, i+1):
            final_rates[j, i] = temp_f[j, i] + temp_res.x[0]

print("*"*100,"\nfinal adjusted rates\n",final_rates, "\n","*"*100)

0.009998219066982433
0.02064955986940258





array([[0.05      , 0.06525676, 0.0817197 ],
       [0.        , 0.05524009, 0.07064956],
       [0.        , 0.        , 0.0615861 ]])

You may also enjoy