# Using Market Prices of Bonds to Model the Discount Curve with a Polynomial (Math + R)

I have a small program I'm building to interpolate the discount curve from a portfolio of benchmark bonds. If anyone has any guesses as to whether it's my process, or my code that's messed up I would greatly appreciate it.

The math behind it is for the most part straight forward compared to most interest rate models:

The market price of a bond is the discounted value of its cash flows:

MPj = Cj1*D1 + Cj2*D2 +… + Cji*Di + Error
Where
MPj = Cash price of bond j
Cji = ith cash flow for bond j
Di = Discount factor for the ith cash flow


Use a polynomial to model the discount factor:

Di = B0 + B1*Ti  + B2*Ti^2 + … + Bk*Ti^k
Where:
Di = Discount factor at time = i
k = The degree of the polynomial
Ti = Time until the ith cash flow, in years
Bk = The coefficients of the model which describe how the time to cash flow determines the discount factor


Substitute the second equation into the first:

MPj = Cj1 [B0 + B1*T1 + B2*T12 + … + Bk*T1k] +… + Cji [B0 + B1*Ti + B2*Ti2 + … + Bk*Tik ] + Error


and in matrix form: http://i.imgur.com/VUwPZ1q.png

From there, take the squared error of this equation, and minimize it by changing the coefficients.

Does this process make sense?

When I try to implement it, I end up with very large squared errors (~73.5 over 15 bonds).

This is the data I'm using: https://www.dropbox.com/s/zeg5xyt5kq4xcpm/Sample%20Data.csv?dl=0

This is the code I'm using. Any comments are appreciated.

library(dplyr)
library(lubridate)

Polynomial_Degree = 3
Start_Date <- as.Date(ymd(19990122))
Face_Value <- 100
Default_Coeffs_Guess <- rep(0,Polynomial_Degree + 1)

Databank <- read.csv(paste(getwd(), "Sample Data.csv", sep="/"), stringsAsFactors = F)

Coupon_Values_List <- Databank %>%
select(COUPON) %>%
data.frame %>%
rename(.,Coupon_Values_List = COUPON)

Coupon_Values_List2 <- Databank %>%
select(COUPON) %>%
transmute(., Coupon_Values_List2 = COUPON +1) %>%
data.frame

Maturity_Date_List <- Databank %>%
select(MATURITY) %>%
lapply(mdy) %>%
lapply(as.Date) %>%
data.frame

Bond_Midpricing <- Databank %>%
as.matrix

#Functions 1:
Coupon_Count_Function <- function(Maturity_Date, Coupon_Payment_Frequency){

Coupon_Count <- 0
Coupon_Count <- as.double(ceiling(((as.Date(Maturity_Date)-as.Date(Start_Date))*Coupon_Payment_Frequency)/365)-1)

return(Coupon_Count)
}

#Function 2:
Coupon_Dates_Function <- function(Maturity_Date, Coupon_Payment_Frequency){

Coupon_Payment_Dates <- NULL
Coupon_Count=Coupon_Count_Function(Maturity_Date,2)

for (i in 0:Coupon_Count)
{x <- as.numeric(((as.Date(as.Date(Maturity_Date)-i*365/Coupon_Payment_Frequency))-Start_Date)/365)
Coupon_Payment_Dates <- rbind(Coupon_Payment_Dates, data.frame(x))}

return((Coupon_Payment_Dates))
}

#Function 3:
Coupon_Time_Func <- function(Maturity_Dates, Coupon_Payment_Frequency) {

All_Coupon_Time_List <- apply(Maturity_Dates, 1, Coupon_Dates_Function, Coupon_Payment_Frequency = 2)

return(All_Coupon_Time_List)
}

#Function 4:
Independent_Variable_Matrix_Func <- function(Maturity_Date_List, Coupon_Payment_Frequency, Coupon_Values_List,
Polynomial_Degree){

df <- data.frame(I(Coupon_Time_Func(Maturity_Date_List, Coupon_Payment_Frequency)),Coupon_Values_List)
names(df) <- c("All_Coupon_Time_List", "Coupon_Values_List")

Temp_List <- vector(mode = "list", Polynomial_Degree+1)
for (i in 1:(Polynomial_Degree+1)) {

Summation_Column_Coupons <- NULL
Summation_Column_Coupons <- mapply(FUN = function(x, y) {y* sum(x^(i-1))}, df$All_Coupon_Time_List, df$Coupon_Values_List)

Summation_Column_With_Principle <- NULL
Summation_Column_With_Principle <- mapply(Summation_Column_Coupons, df\$All_Coupon_Time_List, Face_Value,
FUN = function(x,y,z) {x+max(y^(i-1))*z})

Temp_List[[i]] <- Summation_Column_With_Principle
}

Independent_Variable_Matrix <- NULL
Independent_Variable_Matrix <- as.matrix(t(rbind(Independent_Variable_Matrix, do.call(rbind, Temp_List))))
return(Independent_Variable_Matrix)
}

Independent_Variable_Matrix <- Independent_Variable_Matrix_Func(Maturity_Date_List, Coupon_Payment_Frequency, Coupon_Values_List,
Polynomial_Degree)

#Function 5:
Sum_Error_Squared_Func <- function(Coeffs){

Estimated_Prices_Vector <-  Independent_Variable_Matrix %*% Coeffs

Sum_Error_Squared <- sum((Bond_Midpricing - Estimated_Prices_Vector)^2)

return(Sum_Error_Squared)
}

optim (par = Default_Coeffs_Guess,
fn = Sum_Error_Squared_Func,
gr = NULL,
method  = "L-BFGS-B",
lower   = c(-1000,-1000,-1000,-1000),
upper   = c(1000,1000,1000,1000))