################################################################
#### Example 3.16: Immunization for a non-life insurer 
################################################################

###
# Load data  
###

# Load the Swedish daily zero rates data
data <- read.csv("Data/yieldMatrix20101125.csv", header = FALSE)
data <- data[,2:17]
# Load dates for the Swedish zero rates data
rdates <- read.table("Data/yieldMatrixDates.txt", header = FALSE)
dates <- rdates[,1]
dates <- as.Date(dates,"%d/%m/%y")

# The times to maturity on the zero rates data
# '1W', '1M', '2M', '3M', '6M', '9M', '1Y', '2Y', '3Y', '4Y', '5Y', '6Y', '7Y', '8Y', '9Y', '10Y'

# Extract 3m observations
n <- length(dates)
tdates <- c(dates[1])
tdata <- c(data[1,])
mdates <- seq(dates[1],dates[n],by="3 months")
j <- 2
for(i in 2:n){
   if(months(dates[i])==months(mdates[j]) & 
months(dates[i])!=months(dates[i-1]) & j <= length(mdates))
     {
     tdata <- rbind(tdata,data[i,])
     tdates <- c(tdates,dates[i])
     j <- j+1
   }
}
quadata <- tdata # contains quartely zero rates data
quadates <- tdates # contains corresponding dates

# Problem with missing values on 2003-08-01 and 2003-11-03. solve it by linear interpolation
tmp <- (quadata[29,6]+quadata[29,8])/2
quadata[29,7] <- tmp
tmp <- (quadata[30,6]+quadata[30,8])/2
quadata[30,7] <- tmp

# Problem with missing 10yr observations. Fill in...
quadata[29,16] <- quadata[29,15]+0.073
quadata[30,16] <- quadata[30,15]+0.075

# Define the times to maturity
timetomaturity <- c(0.25,0.5,0.75,1:10)

# Extract zero rates corresponding to the times to maturity
quayielddata <- quadata[,4:16]

# Extract zero rate changes
quayielddatadiff <- diff(as.matrix(quayielddata),lag=1)
n <- length(quayielddatadiff[,1])

###  
# Perform PCA on the zero rate changes
###

dim(quayielddatadiff)
	#[1] 54 13
quayielddatadiff<-quayielddatadiff/100 # Normalize 
quaeigen<- eigen(cov(quayielddatadiff)) 
pc <- quaeigen$vectors #The principal components o1, o2, etc

###
# Define bonds and the liability
###

expliacf<-c(27.20,49.47,67.70,82.63,67.65,55.39,45.35,37.13,30.40,24.89,20.38,16.68,13.66,11.18,9.16,7.50,6.14,5.02,4.11,3.37,2.76,2.26,1.85,1.51,1.24,1.01,0.83,0.68,0.56,0.46,0.37,0.31,0.25,0.20,0.17,0.14,0.11,0.09,0.08,0.06)

# Define function to convert zero rates from the times to maturity to each quarter: use linear interpolation

ydyearstoquarters<-function(ydiffs)
{
tmpvals<-c(1:40)
tmpvals[1:4]<-ydiffs[1:4]

tmpvals[5]<-(3*ydiffs[4]+ydiffs[5])/4
tmpvals[6]<-(2*ydiffs[4]+2*ydiffs[5])/4
tmpvals[7]<-(ydiffs[4]+3*ydiffs[5])/4
tmpvals[8]<-ydiffs[5]

tmpvals[9]<-(3*ydiffs[5]+ydiffs[6])/4
tmpvals[10]<-(2*ydiffs[5]+2*ydiffs[6])/4
tmpvals[11]<-(ydiffs[5]+3*ydiffs[6])/4
tmpvals[12]<-ydiffs[6]

tmpvals[13]<-(3*ydiffs[6]+ydiffs[7])/4
tmpvals[14]<-(2*ydiffs[6]+2*ydiffs[7])/4
tmpvals[15]<-(ydiffs[6]+3*ydiffs[7])/4
tmpvals[16]<-ydiffs[7]

tmpvals[17]<-(3*ydiffs[7]+ydiffs[8])/4
tmpvals[18]<-(2*ydiffs[7]+2*ydiffs[8])/4
tmpvals[19]<-(ydiffs[7]+3*ydiffs[8])/4
tmpvals[20]<-ydiffs[8]

tmpvals[21]<-(3*ydiffs[8]+ydiffs[9])/4
tmpvals[22]<-(2*ydiffs[8]+2*ydiffs[9])/4
tmpvals[23]<-(ydiffs[8]+3*ydiffs[9])/4
tmpvals[24]<-ydiffs[9]

tmpvals[25]<-(3*ydiffs[9]+ydiffs[10])/4
tmpvals[26]<-(2*ydiffs[9]+2*ydiffs[10])/4
tmpvals[27]<-(ydiffs[9]+3*ydiffs[10])/4
tmpvals[28]<-ydiffs[10]

tmpvals[29]<-(3*ydiffs[10]+ydiffs[11])/4
tmpvals[30]<-(2*ydiffs[10]+2*ydiffs[11])/4
tmpvals[31]<-(ydiffs[10]+3*ydiffs[11])/4
tmpvals[32]<-ydiffs[11]

tmpvals[33]<-(3*ydiffs[11]+ydiffs[12])/4
tmpvals[34]<-(2*ydiffs[11]+2*ydiffs[12])/4
tmpvals[35]<-(ydiffs[11]+3*ydiffs[12])/4
tmpvals[36]<-ydiffs[12]

tmpvals[37]<-(3*ydiffs[12]+ydiffs[13])/4
tmpvals[38]<-(2*ydiffs[12]+2*ydiffs[13])/4
tmpvals[39]<-(ydiffs[12]+3*ydiffs[13])/4
tmpvals[40]<-ydiffs[13]

tmpvals
}


# Initial quarterly zero rates and discount factors

initrates<-ydyearstoquarters(as.double(quayielddata[55,]))/100
initdiscvec<-exp(-(1:40)*initrates/4)
#quarterdiscvec<-exp(-((1:40)-1)*initrates/4) 


# Present value of liability
L0<-sum(expliacf*initdiscvec)

# L0 = 582.12

# Cash flow of the four bonds in matrix form

aflowmat<-matrix(0,4,40)
aflowmat[1,1]<-105.25
aflowmat[2,3]<-5.5
aflowmat[2,7]<-105.5
aflowmat[3,2]<-4.5
aflowmat[3,6]<-4.5
aflowmat[3,10]<-4.5
aflowmat[3,14]<-4.5
aflowmat[3,18]<-104.5
aflowmat[4,4]<-5
aflowmat[4,8]<-5
aflowmat[4,12]<-5
aflowmat[4,16]<-5
aflowmat[4,20]<-5
aflowmat[4,24]<-5
aflowmat[4,28]<-5
aflowmat[4,32]<-5
aflowmat[4,36]<-5
aflowmat[4,40]<-105

# To get the four immunization weights we solve a matrix equation

amat<-matrix(0,4,4)
lmat<-matrix(0,4,1)
for(i in (1:4))
{
amat[1,i]<-sum(aflowmat[i,]*initdiscvec)
amat[2,i]<-sum(((1:40))*aflowmat[i,]*initdiscvec*ydyearstoquarters(quaeigen$vectors[,1])/4)
amat[3,i]<-sum(((1:40))*aflowmat[i,]*initdiscvec*ydyearstoquarters(quaeigen$vectors[,2])/4)
amat[4,i]<-sum(((1:40))*aflowmat[i,]*initdiscvec*ydyearstoquarters(quaeigen$vectors[,3])/4)
}
lmat[1,1]<-sum(expliacf*initdiscvec)
lmat[2,1]<-sum(((1:40))*expliacf*initdiscvec*ydyearstoquarters(quaeigen$vectors[,1])/4)
lmat[3,1]<-sum(((1:40))*expliacf*initdiscvec*ydyearstoquarters(quaeigen$vectors[,2])/4)
lmat[4,1]<-sum(((1:40))*expliacf*initdiscvec*ydyearstoquarters(quaeigen$vectors[,3])/4)


imweights<-solve(amat)%*%lmat # The optimal portfolio weights 
imweights

# > imweights
#           [,1]
#[1,] 1.11819503
#[2,] 3.74204295
#[3,] 0.46077515
#[4,] 0.07171441

###
# Simulation study
###


library(MASS)

normsampledata<-mvrnorm(100000,colMeans(quayielddatadiff,dim=1),cov(quayielddatadiff))
L1vals<-(1:100000)
A1vals<-(1:100000)
for(i in (1:100000))
{
L1vals[i]<-sum(expliacf*exp(-((1:40)-1)*(initrates+ydyearstoquarters(normsampledata[i,]))/4))
A1vals[i]<-sum((imweights[1]*aflowmat[1,]+imweights[2]*aflowmat[2,]+imweights[3]*aflowmat[3,]+imweights[4]*aflowmat[4,])*exp(-((1:40)-1)*(initrates+ydyearstoquarters(normsampledata[i,]))/4))
}


# Histograms of hedging errors
# Hedging error without immunization
hist(L0/initdiscvec[1]-L1vals,200,xlab="",ylab="",main="",yaxt="n",col="black")

# Hedging error with immunization
hist(A1vals-L1vals,200,xlab="",ylab="",main="",yaxt="n",col="black")

# Predicted cash flow from liability (negative) and 
# cash flow from the mathching bond portfolio (positive)

plot((1:40)/4,imweights[1]*aflowmat[1,1:40]+imweights[2]*aflowmat[2,1:40]+imweights[3]*aflowmat[3,1:40]+imweights[4]*aflowmat[4,1:40],type="h",xlab="",ylab="",ylim=c(-max(expliacf),max(imweights[2]*aflowmat[2,1:40])))
par(new=T)
plot((1:40)/4,-expliacf[1:40],type="h",xlab="",ylab="",ylim=c(-max(expliacf),max(imweights[2]*aflowmat[2,1:40])))
lines(c(0,10),c(0,0))