###############################################################
### Example 7.12: Sampling from the predictive distribution (chain ladder)
###############################################################

# Interest rates

# 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(1:10)

# Extract zero rates corresponding to the times to maturity
quayielddata <- quadata[,7:16]
zerorates <- as.double(quayielddata[55,])

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

###
# Define insurance parameters
###

# The upper triangle of claim amounts

claimtri <- rbind(c(357848, 766940, 610542, 482940, 527326, 574398, 146342, 139950, 227229, 67948), c(352118, 884021, 933894, 1183289, 445745, 320996, 527804, 266172, 425046,0), c(290507, 1001799, 926219, 1016654, 750816, 146923, 495992, 280405,0,0), c(310608, 1108250, 776189, 1562400, 272482, 352053, 206286,0,0,0), c(443160, 693190, 991983, 769488, 504851, 470639,0,0,0,0), c(396132, 937085, 847498, 805037, 705960,0,0,0,0,0), c(440832, 847631, 1131398, 1063269,0,0,0,0,0,0), c(359480, 1061648, 1443370,0,0,0,0,0,0,0), c(376686, 986608,0,0,0,0,0,0,0,0), c(344014,0,0,0,0,0,0,0,0,0))

n <- dim(claimtri)[1]
m <- dim(claimtri)[2]


###
# Extracting residuals
###

# Convert the claim amounts to cumulative amounts
cumtri <- matrix(0,m,n)
for(i in 1:m){
	cumtri[i,1:(n-i+1)] <- cumsum(claimtri[i,1:(n-i+1)])
	}
cumtri

# Compute the development factors 
dev.fact <- c()
for(j in 1:(n-1)){
	dev.fact <- c(dev.fact,sum(cumtri[1:(m-j),(j+1)])/sum(cumtri[1:(m-j),j]))
	}	
dev.fact <- c(dev.fact,1)
dev.fact

# Compute estimates of the mean of claim amounts
newcumtri <- matrix(0,m,n)
for(i in 1:(m-1)){
	newcumtri[i,n-i+1] <- cumtri[i,n-i+1]
	for(j in (n-i):1){
		newcumtri[i,j] <- newcumtri[i,(j+1)]/dev.fact[j]
		}
}
newcumtri[m,1] <- cumtri[m,1]
fitclaimtri <- newcumtri
for(i in 1:(m-1)){
		fitclaimtri[i,2:(n-i+1)] <- diff(newcumtri[i,1:(n-i+1)])
	}


# Compute residuals
restri <- (claimtri - fitclaimtri)/fitclaimtri


# Store residuals in a vector
resvec <- c()
for(i in 1:m){
	resvec <- c(resvec,restri[i,1:(n-i+1)])
	}



###
# Compute time 0  value of assets and liabilities
###

# Compute predictions of the cumulative lower triangle
cumclaimpred <- cumtri
	for(i in 2:m){
		cumclaimpred[i,(n-i+2):n] <- cumtri[i,(n-i+1)]*cumprod(dev.fact[(n-i+1):(n-1)]) 
		}
cumclaimpred <- rbind(cumclaimpred,mean(fitclaimtri[,1])*c(1,cumprod(dev.fact[1:n-1])))

claimpredlow <- cumclaimpred
for(i in 1:m+1){
	claimpredlow[i,2:n] <- diff(cumclaimpred[i,])
}


payvec0 <- 0*c(1:n)
	for(i in 2:m+1){
		for(j in (n-i+2):n){
			payvec0[i+j-11] <- payvec0[i+j-11] + claimpredlow[i,j] 
			}
 	}
 	payvec0 #predicted payment amounts for each year
 	
 L0 <- payvec0%*%exp(-as.double(zerorates)*c(1:10)/100) # Value of liability
 A0 <- L0 # value of asset portfolio with weights payvec0

###
# Sampling the predictive distribution
###

lowtrizero <- matrix(0,n,m)
for(i in 1:n){
	for(j in 1:(n-i+1)){
		lowtrizero[i,j] <- 1
		}
	}


N <- 10000
paymat <- matrix(0,N,10)

L1vec <- c()
A1vec <- c()

neg.rate <- 0

for(k in 1:N){
	m <- dim(fitclaimtri)[1]
	n <- dim(fitclaimtri)[2]
	
	restri.b <- sample(resvec, n*m,replace = TRUE)
	dim(restri.b) <- c(m,n)
	restri.b <- restri.b*lowtrizero
	claimtri.b <- fitclaimtri*(1+restri.b) # Sampled upper tri
	
	cumtri.b <- matrix(0,m,n) # Sampled cumulative tri
	for(i in 1:m){ 
		cumtri.b[i,1:(n-i+1)] <- cumsum(claimtri.b[i,1:(n-i+1)])
	}
	
	dev.fact.b <- c() # Development factors of the sampled tri
	for(j in 1:(n-1)){
		dev.fact.b <- c(dev.fact.b,sum(cumtri.b[1:(m-j),(j+1)])/sum(cumtri.b[1:(m-j),j]))
	}	
	dev.fact.b <- c(dev.fact.b,1)
	
	cumtri.bb <- rbind(cumtri.b,0*c(1:n)) # Sample new diagonal
	for(i in 2:m){
		cumtri.bb[i,(n-i+2)] <- cumtri.b[i,(n-i+1)]*(dev.fact.b[n-i+1]-1)*(1+sample(resvec,1,replace=TRUE)) + cumtri.b[i,(n-i+1)]
		}
		cumtri.bb[n+1,1] <- mean(cumtri.b[,1])*(1+sample(resvec,1))
		
	dev.fact.bb <- c() # Update the development factors
	for(j in 1:(n-1)){
		dev.fact.bb <- c(dev.fact.bb,sum(cumtri.bb[1:(m+1-j),(j+1)])/sum(cumtri.bb[1:(m+1-j),j]))
	}	
	dev.fact.bb <- c(dev.fact.bb,1)
	cumpredlow.bb <- cumtri.bb
	for(i in 3:(m+1)){
		cumpredlow.bb[i,(n-i+3):n] <- cumpredlow.bb[i,(n-i+2)]*cumprod(dev.fact.bb[(n-i+2):(n-1)]) 
		}
	claimpredlow.bb <- cumpredlow.bb
	for(i in 1:(m+1)){
		claimpredlow.bb[i,2:n] <- diff(cumpredlow.bb[i,])
		}
		
	payvec.bb <- 0*c(1:n)
	for(i in 2:(m+1)){
		for(j in (n-i+2):n){
			payvec.bb[i+j-11] <- payvec.bb[i+j-11] + claimpredlow.bb[i,j] 
			}
 	}	
 	paymat[k,] <- payvec.bb
 	zrindex <- sample(c(1:qn),4, replace = TRUE)
 	zeroratesscenario <- (zerorates + colSums(quayielddatadiff[zrindex,]))/100
 	#zeroratesscenario <- (zerorates + quayielddatadiff[zrindex,])/100
 	zerorates.b <- as.double(c(0,zeroratesscenario[1:9]))
 	L1vec <- c(L1vec,paymat[k,]%*%exp(-zerorates.b*c(0:9)))	
 	A1vec <- c(A1vec, payvec0%*%exp(-zerorates.b*c(0:9)))
 	if(min(zerorates.b) < 0) neg.rate <- neg.rate+1
	}

# Note: with the procedure outlined above interest rates generated by historical simulation sometimes become negative!!! The variable 'neg.rate' counts how many simulations that contain negative interest rates. The reason is that historically interest rates has been declining and is at a very low level. A better alternative is to use a different distribution for future interest rates, possibly adding a positive mean, or centering around the forward rates.   


hist(L1vec,breaks = 100, col="black")
points(L0*exp(zerorates[1]/100),0,col="red", pch = 19)

hist(A1vec, breaks = 75, col="black")
points(A0*exp(zerorates[1]/100),0,col="red", pch=19)

hist(A1vec-L1vec,breaks = 75, col="black",xlab="",ylab="", axes = FALSE, main="")
axis(1)