R function for computing the BDT short rate tree
# Info: This function is based on the code by Zhang Lei in "A Java Applet for pricing Bonds and Bond Options using the Black-Derman-Toy model", master thesis, 2005.
# Description: Constructs the short rate tree using the data calibration model proposed by Black, Derman and Toy.
# Parameter 1 - yields: market zero-coupon bond yields.
# Parameter 2 - volatilities: volatilities in force for each of the zero-coupon bonds at time 1. Note that volatilites[1] is undefined.
# Parameter 3 - dt: time period of each step expressed in years.
# Example 1: BDT(yields = c(0.10, 0.11, 0.12, 0.125, 0.13), volatilities = c(NA, 0.19, 0.18, 0.17, 0.16), dt = 1)
# Example 2: BDT(yields=c(0.2004800, 0.1924200, 0.1891550, 0.1858900, 0.1836150, 0.1813400, 0.1790650, 0.1767900, 0.1741775, 0.1715650, 0.1689525, 0.1663400, 0.1644950, 0.1626500, 0.1608050, 0.1589600, 0.1586400, 0.1583200, 0.1580000, 0.1576800), volatilities=c(NA, 0.0393090, 0.0374260, 0.0370270, 0.0358750, 0.0347230, 0.0335710, 0.0324190, 0.0324168, 0.0324145, 0.0324123, 0.0324100, 0.0324275, 0.0324450, 0.0324625, 0.0324800, 0.0318398, 0.0311995, 0.0305593, 0.0299190), dt=0.25)
BDT <- function(yields, volatilities, dt) {
N = length(yields)
sqdt = sqrt(dt)
r <- matrix(numeric(0), (2*N)-1, N)
r[N,1] = yields[1]
d <- matrix(numeric(0), (2*N)-1, N)
d[N,1] = (1+yields[1])^(-dt)
# Calculate price of pure discount bonds
P <- numeric(0)
for(i in 1:N) {
P[i] = (1+yields[i])^(-i*dt)
}
Pu <- numeric(0)
Pu[1] = 1
Pd <- numeric(0)
Pd[1] = 1
epsilon = 1e-07
expfunc <- numeric(0)
# Calculate Pu and Pd
for (i in 2:N) {
# Calculate Pu with Newton-Raphson
Pu[i] = Pu[i-1]
expfunc[i] = exp(-2*volatilities[i]*sqdt)
tempFunc = Pu[i]+(1 - expfunc[i] + expfunc[i]*Pu[i]^(-1/((i-1)*dt)))^(-(i-1)*dt) - 2*P[i]*(1+r[N,1])^dt
while (abs(tempFunc) > epsilon) {
tempFuncDer = 1 + (1 - expfunc[i] + expfunc[i]*Pu[i]^(-1/((i-1)*dt)))^(-(i-1)*dt-1)*expfunc[i]*Pu[i]^(-1/((i-1)*dt)-1)
Pu[i] = Pu[i] - tempFunc/tempFuncDer
tempFunc = Pu[i]+(1 - expfunc[i] + expfunc[i]*Pu[i]^(-1/((i-1)*dt)))^(-(i-1)*dt) - 2*P[i]*(1+r[N,1])^dt
}
# Calculate Pd with formula
Pd[i] = (1 - expfunc[i] + expfunc[i]*Pu[i]^(-1/((i-1)*dt)))^(-(i-1)*dt)
}
U <- numeric(0)
sig <- numeric(0)
for (i in 2:N) {
U[i-1] = yields[i]
sig[i-1] = volatilities[i]
}
Qu <- matrix(numeric(0), (2*N)-1, N-1)
Qu[N+1,1] = 1
Qd <- matrix(numeric(0), (2*N)-1, N-1)
Qd[N-1,1] = 1
# Build short rate tree
for (i in 1:(N-1)) {
# Build Qu and Qd trees
if (i > 1) {
Qu[N+i,i] = 0.5*Qu[N+i-1,i-1]*d[N+i-1,i] # top row of Qu
Qu[N+2-i,i] = 0.5*Qu[N+2-i+1,i-1]*d[N+2-i+1,i] # low row of Qu
Qd[N-2+i,i] = 0.5*Qd[N-2+i-1,i-1]*d[N-2+i-1,i] # top row of Qd
Qd[N-i,i] = 0.5*Qd[N-i+1,i-1]*d[N-i+1,i] # low row of Qd
if (i > 2) {
for (j in seq(N+2-i+2, N+i-2, 2)) { # inner rows of Qu
Qu[j,i] = 0.5*Qu[j-1,i-1]*d[j-1,i] + 0.5*Qu[j+1,i-1]*d[j+1,i]
}
for (j in seq(N-i+2, N-2+i-2, 2)) { # inner rows of Qd
Qd[j,i] = 0.5*Qd[j-1,i-1]*d[j-1,i] + 0.5*Qd[j+1,i-1]*d[j+1,i]
}
}
}
# Calculate U and sig using Newton-Raphson
# calculate initial values
expfunc <- numeric(0)
for (j in seq(N-i, N+i, 2)) {
expfunc[j] = exp(sig[i]*(j-N)*sqdt)
}
sumtemp = 0
for (j in seq(N-i+2, N+i, 2)) {
sumtemp = sumtemp + Qu[j,i]/(1+U[i]*expfunc[j])^dt
}
tempFunc = sumtemp - Pu[i+1]
sumtemp = 0
for (j in seq(N-i, N+i-2, 2)) {
sumtemp = sumtemp + Qd[j,i]/(1+U[i]*expfunc[j])^dt
}
tempFunc2 = sumtemp - Pd[i+1]
while ( (abs(tempFunc) > epsilon) | (abs(tempFunc2) > epsilon)) {
# calculate derivatives
sumtempDer = 0
for (j in seq(N-i+2, N+i, 2)) {
sumtempDer = sumtempDer - Qu[j,i]*expfunc[j]*dt/((1+U[i]*expfunc[j])^(dt+1))
}
tempFuncDer = sumtempDer
sumtempDer = 0
for (j in seq(N-i+2, N+i, 2)) {
sumtempDer = sumtempDer - Qu[j,i]*U[i]*dt*(j-N)*sqdt*expfunc[j]/((1+U[i]*expfunc[j])^(dt+1))
}
tempFuncDer2 = sumtempDer
sumtempDer = 0
for (j in seq(N-i, N+i-2, 2)) {
sumtempDer = sumtempDer - Qd[j,i]*expfunc[j]*dt/((1+U[i]*expfunc[j])^(dt+1))
}
tempFunc2Der = sumtempDer
sumtempDer = 0
for (j in seq(N-i, N+i-2, 2)) {
sumtempDer = sumtempDer - Qd[j,i]*U[i]*dt*(j-N)*sqdt*expfunc[j]/((1+U[i]*expfunc[j])^(dt+1))
}
tempFunc2Der2 = sumtempDer
# calculate new values for U and sig
jacob = (tempFuncDer*tempFunc2Der2) - (tempFunc2Der*tempFuncDer2)
U[i] = U[i] - ((tempFunc*tempFunc2Der2 - tempFunc2*tempFuncDer2)/jacob)
sig[i] = sig[i] - ((tempFunc2*tempFuncDer - tempFunc*tempFunc2Der)/jacob)
# update functions
for (j in seq(N-i, N+i, 2)) {
expfunc[j] = exp(sig[i]*(j-N)*sqdt)
}
sumtemp = 0
for (j in seq(N-i+2, N+i, 2)) {
sumtemp = sumtemp + Qu[j,i]/(1+U[i]*expfunc[j])^dt
}
tempFunc = sumtemp - Pu[i+1]
sumtemp = 0
for (j in seq(N-i, N+i-2, 2)) {
sumtemp = sumtemp + Qd[j,i]/(1+U[i]*expfunc[j])^dt
}
tempFunc2 = sumtemp - Pd[i+1]
}
# Compute short rates
for (j in seq(N-i, N+i, 2)) {
r[j,i+1] = U[i]*exp(sig[i]*(j-N)*sqdt)
d[j,i+1] = (1+r[j,i+1])^(-dt)
}
}
r
}