library(tseries)
library(lubridate)

### In this example, we will use the following 5 bonds ###

Bond = c("Gov Bond 2018",
         "Gov Bond 2020",
         "Gov Bond 2023",
         "Gov Bond 2027",
         "Gov Bond 2039")
Bond

FaceValue = rep(100, length(Bond))

Coupon = c(0.25, 0.25, 1.50, 0.50, 4.50) #bond coupons

cleanPrice = c(100.960,102.285,109.700,99.993,168.794) 

Maturity = as.Date(c("2018-11-15","2020-11-15", "2023-11-15","2027-11-15","2039-11-15"))

Maturity
bondData <- data.frame(Bond, FaceValue, Coupon, cleanPrice, Maturity)
bondData
couponsToMaturity <- year(bondData$Maturity)-year("2017-10-16") + 1 ## + 1 because we want to have the coupon payment for 2017 in the cashflows

cashFlowDates <- c(as.Date("2017-10-16"), seq.Date(from= as.Date("2017-11-15"), to = as.Date("2039-11-15"), by = 'year'))
cashFlowDates

# Calculate Dirty Price (Weights, if budget is 100000)
bondData$dirtyPrice <- bondData$cleanPrice + bondData$Coupon * as.numeric(difftime("2017-10-16", "2016-11-15", units = 'days'))/365

bondData$NumberofBonds = rep(100000, length(Bond))/bondData$dirtyPrice
#Ç100000/988.2456
# CashFlow for Dirty Price
# Create a CashFlow Matrix
cashFlowsDirty <- matrix(0L, nrow = length(Bond), ncol = length(cashFlowDates))
rownames(cashFlowsDirty) <- Bond
colnames(cashFlowsDirty) <- as.character(cashFlowDates)
cashFlowsDirty

for (i in 1:length(Bond)){
  cashFlowsDirty[i,1] <- -bondData$dirtyPrice[i]   # Negative Cashflow when bond is bought
  for (j in 1:couponsToMaturity[i]){
    cashFlowsDirty[i,j+1] <- bondData$Coupon[i]   # Coupon payment
    if (j == couponsToMaturity[i]) {
      cashFlowsDirty[i,j+1] <- cashFlowsDirty[i,j+1] + bondData$FaceValue[i]  # Face value of bond paid at Maturity
    }
  }
}

cashFlowsDirty
bondData$NumberofBonds
cashFlowsDirty

cashFlowsDirtyPortfolio = t(bondData$NumberofBonds) %*% cashFlowsDirty
cashFlowsDirtyPortfolio
# Create bond valuation function
# We need to take into account that it is not a year until the first coupon paypent,
bval <- function(i, cf, t=seq(along = cf)-1){
  t[2:length(t)] <- t[2:length(t)] - ( 1 - as.numeric(difftime("2017-11-15", "2017-10-16", units='days'))/365)  # here we account for the time until coupon payment
  sum(cf / (1 + i)^t)
}


# Calculate Yield to Maturity
YTMdirty = c()
for (i in 1:length(Bond)){
  YTMdirty[i] <- uniroot(bval, interval = c(-0.5, 2), cf = cashFlowsDirty[i,])$root
}
YTMdirty
bondData$YTMdirty <- YTMdirty*100
bondData

YTMdirtyportfolio <- uniroot(bval, interval = c(-0.5, 2), cf = cashFlowsDirtyPortfolio)$root
YTMdirtyportfolio <- YTMdirtyportfolio*100
YTMdirtyportfolio
# CashFlow for Dirty Price
CashFlowClean <- cashFlowsDirty
CashFlowClean[,1] <- -bondData$cleanPrice
CashFlowClean[,2] <- CashFlowClean[,2] - bondData$Coupon *  (365-as.numeric(difftime("2017-11-15", "2017-10-16", units = 'days')))/365

CashFlowClean

# Calculate Yield to Maturity
YTMclean = c()
for (i in 1:length(Bond)){
  YTMclean[i] <- uniroot(bval, interval = c(-0.5, 2), cf = CashFlowClean[i,])$root
}

bondData$YTMclean <- YTMclean*100
bondData


# Function for Calculating Duration
calcDuration <- function(Price, cf, t, r){
  1/Price * sum( t * cf / (1 + r)^(t+1))
}

# Function For Calculating Convexity
calcConvexity <- function(Price, cf, t, r){
  1/Price * sum( t * (t+1) * cf / (1 + r)^(t+2))
}

dur <- c() # Empty array to store Duration Calculation
conv <- c() # Empty array to store Convexity Calculation

for (i in 1:length(bondData$Bond)){
  p <- bondData$dirtyPrice[i]
  cf <- cashFlowsDirty[i, ]
  
  t <- seq(along = cf)-1
  t[2:length(t)] <- t[2:length(t)] - ( 1 - as.numeric(difftime("2017-11-15", "2017-10-16", units='days'))/365)  # here we account for the time until coupon payment
  
  r <- bondData$YTMdirty[i]/ 100 
  
  dur <- c(dur, calcDuration(p, cf, t, r))
  conv <- c(conv, calcConvexity(p, cf, t, r))
}

bondData$Duration <- dur
bondData$Convexity <- conv

bondData

#Note that amount invested is a constant at 100  per bond
# Portfolio Duration Calculation

sum(bondData$Duration)/5
DurationPortfolio <- calcDuration(500000, cashFlowsDirtyPortfolio, t , YTMdirtyportfolio/100)
DurationPortfolio

sum(bondData$Convexity)/5
ConvexityPortfolio <- calcConvexity(500000, cashFlowsDirtyPortfolio, t, YTMdirtyportfolio/100)
ConvexityPortfolio

### We now compute the potential market decline by Taylor Series Approximation ###
# A 120 basis point change is equivalent to a change by 1.2%
# Change in market value = duration x change in rate + convexity x (change in rate)^2 / 2
valuechange <- (-DurationPortfolio*0.012 + ConvexityPortfolio*(0.012)^2 /2)
valuechange

PriceBP=sum(bondData$dirtyPrice*bondData$NumberofBonds)
PriceBP

delta_PBP=valuechange*PriceBP
delta_PBP
New_PriceBP=PriceBP+delta_PBP
New_PriceBP