### Non-Life Insurance: Mathematics and Statistics ### Exercise Sheet 8 ### Exercise 1 ### Define the function KK_premium with the variables: ### lambda = mean number of claims ### mu = mean parameter of log-normal distribution ### sigma2 = variance parameter of log-normal distribution ### span = span size used in the Panjer algorithm ### shift = shift of the translated log-normal distribution KK_premium <- function(lambda, mu, sigma2, span, shift){ ### We will calculate the distribution of S until M (M = 2500 + 7000) M <- 9500 ### Number of steps m <- floor(M/span) ### We won't have any mass until we reach shift, which happens at the k0-th step k0 <- shift/span ### Initialize array where mass is put to the lower end of the interval g_min <- array(0, dim=c(m+1,1)) ### Initialize array where mass is put to the upper end of the interval g_max <- array(0, dim=c(m+1,1)) ### Discretize the log-normal distribution putting the mass to the lower end of the interval for (k in (k0+1):(m+1)){g_min[k,1] <- pnorm(log((k-k0)*span), mean=mu, sd=sqrt(sigma2))-pnorm(log((k-k0-1)*span), mean=mu, sd=sqrt(sigma2))} ### Discretize the log-normal distribution putting the mass to the upper end of the interval g_max[2:(m+1),1] <- g_min[1:m,1] ### Initialize matrix, where we will store the probability distribution of S f1 <- matrix(0, nrow=m+1, ncol=3) ### Store the probability of getting zero claims (in both lower bound and upper bound) f1[1,1] <- exp(-lambda*(1-g_min[1,1])) f1[1,2] <- exp(-lambda*(1-g_max[1,1])) ### Calculate the values "l * g_{l}" of the discretized claim sizes (lower bound and upper bound), we need these values in the Panjer algorithm h1 <- matrix(0, nrow=m, ncol=3) for (i in 1:m){ h1[i,1] <- g_min[i+1,1]*(i+1) h1[i,2] <- g_max[i+1,1]*(i+1) } ### Panjer algorithm (note that in the Poisson case we have a = 0 and b = lambda*v, which is just lambda here) for (r in 1:m){ f1[r+1,1] <- lambda/r*(t(f1[1:r,1])%*%h1[r:1,1]) f1[r+1,2] <- lambda/r*(t(f1[1:r,2])%*%h1[r:1,2]) f1[r+1,3] <- r * span } ### Maximal and minimal franchise m1 <- 2500 m0 <- 300 ### Number of iterations needed to get to m1 and m0 i1 <- floor(m1/span+1) i0 <- floor(m0/span+1) ### Calculate the part that the insured pays by himself franchise <- array(NA, c(i1, 3)) for (i in i0:i1){ franchise[i,1] <- f1[i,3] ### this represents the franchise franchise[i,2] <- sum(f1[1:i,1]*f1[1:i,3]) + f1[i,3] * (1-sum(f1[1:i,1])) franchise[i,2] <- franchise[i,2] + sum(f1[(i+1):floor(i+7000/span),1]*f1[2:floor(7000/span+1),3])*0.1 + 700 * (1-sum(f1[1:floor(i+7000/span),1])) franchise[i,3] <- sum(f1[1:i,2]*f1[1:i,3]) + f1[i,3] * (1-sum(f1[1:i,2])) franchise[i,3] <- franchise[i,3] + sum(f1[(i+1):floor(i+7000/span),2]*f1[2:floor(7000/span+1),3])*0.1 + 700 * (1-sum(f1[1:floor(i+7000/span),2])) } ### Calculate the price of the monthly premium price <- array(NA, c(i1, 3)) price[,1] <- franchise[,1] ### this represents the franchise price[,2:3] <- (lambda*(exp(mu+sigma2/2)+shift) - franchise[,2:3])/12 price } ### Load the packages stats and MASS require(stats) require(MASS) ### Determine values for the input parameters of the function KK_premium lambda <- 1 mu <- 7.8 sigma2 <- 1 span <- 10 shift <- 100 ### Run the function KK_premium price <- KK_premium(lambda, mu, sigma2, span, shift) ### Plot the monthly pure risk premium as a function of the franchise plot(x=price[,1], y=price[,2], lwd=2, col="blue", type='l', ylab="Monthly pure risk premium", xlab="Franchise", main="Monthly pure risk premium", cex.lab=1.25, cex.main=1.25, cex.axis=1.25) points(x=c(300,500, 1000, 1500, 2000, 2500), y=price[c(300,500, 1000, 1500, 2000, 2500)/span+1,3], pch=19, col="orange") lines(x=c(300,300), y=c(0,price[300/span+1,3]),lty=3, lwd=1.5, col="darkgray") lines(x=c(500,500), y=c(0,price[500/span+1,3]),lty=3, lwd=1.5, col="darkgray") lines(x=c(1000,1000), y=c(0,price[1000/span+1,3]),lty=3, lwd=1.5, col="darkgray") lines(x=c(1500,1500), y=c(0,price[1500/span+1,3]),lty=3, lwd=1.5, col="darkgray") lines(x=c(2000,2000), y=c(0,price[2000/span+1,3]),lty=3, lwd=1.5, col="darkgray") lines(x=c(2500,2500), y=c(0,price[2500/span+1,3]),lty=3, lwd=1.5, col="darkgray") ### Give the monthly pure risk premiums for the six franchises listed on the exercise sheet round(price[floor(c(300, 500, 1000, 1500, 2000, 2500)/span+1),2]) round(price[floor(c(300, 500, 1000, 1500, 2000, 2500)/span+1),3])