#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

save_plots = FALSE

set.seed(0) 

# P 1:
#
bondvalue = function(c, T, r, par){
    #
    # Computes bv = bond values (current prices) corresponding
    # to all values of yield to maturity in the input vector r
    #
    # INPUT
    # c = coupon payment (semiannual)
    # T = time to maturity (in years)
    # r = vector of yields to maturity (semiannual rates)
    # par = par value
    #
    bv = c/r + (par - c/r) * (1+r)^(-2*T)
    bv 
}

# Computes the yield to maturity of a bond paying semiannual coupon payments
#
# price, coupon payment, and time to maturity (in years) are set below
#
# Uses the function "bondvalue"
#
price = 1200 # current price of the bond
C = 40       # coupon payment
T = 30       # time to maturity
par = 1000   # par value of the bond

r = seq( 0.02, 0.05, length.out=500 )
value = bondvalue(C, T, r, par)
yield2M = spline(value, r, xout=price) # spline interpolation

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter3/chap_3_prob_1_plot.eps", onefile=FALSE, horizontal=FALSE) }
plot( r, value, xlab='yield to maturity', ylab='price of bond', type='l', main='par= 1000, coupon payment= 40, T= 30', lwd=2 )
abline(h=price)
abline(v=yield2M$y)
if( save_plots ){ dev.off() }


# P 2:
#
uniroot( function(r) r^2 - 0.5, c(0.7,0.8) )

# P 3:
#
uniroot( function(r) bondvalue(C, T, r, par) - price, c(0.03,0.04) )

# P 4:
#
uniroot( function(r) bondvalue(280, 8, r, 10000) - 9800, c(0.01,0.04) )

# P 5:
#
uniroot( function(r) bondvalue(35, 20, r, 1000) - 1050.0, c(0.01, 0.04) )

# P 6:
#
uniroot( function(c) bondvalue(c, 5, 0.035, 1000) - 950.10, c(10, 40) )

# P 7:
#
mk.zero2 = read.csv('../../BookCode/Data/mk.zero2.csv', header=TRUE)
mk.maturity = read.csv('../../BookCode/Data/mk.maturity.csv', header=TRUE)

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter3/chap_3_prob_7_short_end.eps", onefile=FALSE, horizontal=FALSE) }
plot(mk.maturity[,1], mk.zero2[5,2:56], type="l", xlim=c(0,3), xlab="maturity (monthly)", ylab="yield", main="short end" )
lines(mk.maturity[,1], mk.zero2[6,2:56], lty=2, type="l")
lines(mk.maturity[,1], mk.zero2[7,2:56], lty=3, type="l")
lines(mk.maturity[,1], mk.zero2[8,2:56], lty=4, type="l")
legend("topleft", c("1985-12-01", "1986-01-01", "1986-02-01", "1986-03-01"), lty=1:4)
if( save_plots ){ dev.off() }


if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter3/chap_3_prob_7_long_end.eps", onefile=FALSE, horizontal=FALSE) }
plot(mk.maturity[,1], mk.zero2[5,2:56], type="l", xlim=c(20,30), xlab="maturity (monthly)", ylab="yield", main="long end" )
lines(mk.maturity[,1], mk.zero2[6,2:56], lty=2, type="l")
lines(mk.maturity[,1], mk.zero2[7,2:56], lty=3, type="l")
lines(mk.maturity[,1], mk.zero2[8,2:56], lty=4, type="l")
legend("bottomleft", c("1985-12-01", "1986-01-01", "1986-02-01", "1986-03-01"), lty=1:4)
if( save_plots ){ dev.off() }


# P 8:
#
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter3/chap_3_prob_8.eps", onefile=FALSE, horizontal=FALSE) }

indx = which(mk.zero2[, 1] == '1986-12-01')
plot(mk.maturity[,1], mk.zero2[indx,2:56], type="l", xlim=c(-1, 30), xlab="maturity (monthly)", ylab="yield")

indx = which(mk.zero2[, 1] == '1987-01-01')
lines(mk.maturity[,1], mk.zero2[indx,2:56], type="l", lty=2)

indx = which(mk.zero2[, 1] == '1987-02-01')
lines(mk.maturity[,1], mk.zero2[indx,2:56], type="l", lty=3)

indx = which(mk.zero2[, 1] == '1987-03-01')
lines(mk.maturity[,1], mk.zero2[indx,2:56], type="l", lty=4)

grid()

legend("bottomright", c("1986-12-01", "1987-01-01", "1987-02-01", "1987-03-01"), lty=1:4)

if( save_plots ){ dev.off() }


# P 9:
#
get_forward_rate = function(indx){
    intForward = mk.maturity[, 1] * mk.zero2[indx, 2:56] # T y_T
    xout = seq(0, 20, length=200)
    z1 = spline(mk.maturity[, 1], intForward, xout=xout)
    forward = diff(z1$y) / diff(z1$x) # r(t)
    T_grid = (xout[-1] + xout[-200])/2

    return(list(forward=forward, T_grid=T_grid))
}

indx = 6
fr = get_forward_rate(indx)
plot(fr$T_grid, fr$forward, type='l', lwd=2, ylim=c(0.06, 0.11))
grid()


# Do the problem:
# 
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter3/chap_3_prob_9.eps", onefile=FALSE, horizontal=FALSE) }

indx = which(mk.zero2[, 1] == '1985-12-01')
fr = get_forward_rate(indx)
plot(fr$T_grid, fr$forward, type='l', ylim=c(0.06, 0.115), xlab='maturity (month)', ylab='forward rate')

indx = which(mk.zero2[, 1] == '1986-01-01')
fr = get_forward_rate(indx)
lines(fr$T_grid, fr$forward, type='l')

indx = which(mk.zero2[, 1] == '1986-02-01')
fr = get_forward_rate(indx)
lines(fr$T_grid, fr$forward, type='l')

indx = which(mk.zero2[, 1] == '1986-03-01')
fr = get_forward_rate(indx)
lines(fr$T_grid, fr$forward, type='l')

grid()

legend("bottomright", c("1986-12-01", "1987-01-01", "1987-02-01", "1987-03-01"), lty=1:4)

if( save_plots ){ dev.off() }