demo/MultiChainLadder.R

In mages/ChainLadder: Statistical Methods and Models for Claims Reserving in General Insurance

## MultiChainLadder demos
## Author: Wayne (Yanwei) Zhang, March 2010

# Reproduce results in Zhang (2010)
# Data auto consists of three triangles, paid personal auto, incurred personal auto
# and paid commercial auto. The paid and incurred triangles from personal auto will
# result in divergent Paid-to-incurred ratios under SCL. MCL is inadequate since the
# residual plots have clear downward pattern. GMCL provides a nice solution here to 
# resolve several problems with SCL and MCL. The following illustrates the use of 
# different models step by step.   

auto <- as(auto, "triangles")

# split data so that SCL is used for years 7:10 in MCL and GMCL
da1 <- auto[, 1:7]
da2 <- auto[, 7:10]

# MCL actually fails without using SCL in years 8:9

fit.mcl8 <- MultiChainLadder(auto[, 8:9], model = "MCL", extrap = FALSE)
coef(fit.mcl8)
residCov(fit.mcl8)	

# This is because the input residual covariance is almost singular
if (FALSE){
	inputCov <- fit.mcl8$models[[1]]$residCovEst
	solve(inputCov)
}

#1. fit with SCL 
fit.scl <- MultiChainLadder(auto, "OLS")

## 2. fit with MCL 
fit.mcl <- MultiChainLadder2(auto)

# 3. GMCL with only intercepts
fit.int <- MultiChainLadder2(auto, type = "MCL+int")

### 4. GMCL fit with full models

coefr <- matrix(0, 12, 8)
pos <- cbind(c(1:3, 5:7, 9, 12), 1:8)
coefr[pos] <- 1
int <- 1:6
restrict.regMat <- c(rep(list(coefr), 6),rep(list(NULL), 3))

fit1 <- MultiChainLadder(da1, "SUR", int = int,
			restrict.regMat = restrict.regMat, model = "GMCL")
fit2 <- MultiChainLadder(da2, "OLS")
fit <- Join2Fits(fit1, fit2)
pred <- predict(fit)
mse <- Mse(fit, pred)
fit.full <- JoinFitMse(fit, mse)

# 5. GMCL with parameter selection
# we leave out the detail of the stepwise process 
# and only show the final selected model

# specify differnt restriction matrix for each period
coefr1 <- matrix(0, 12, 6)
pos1 <- cbind(c(2, 3, 5, 7, 9, 12), 1:6)
coefr2 <- matrix(0, 12, 5)
pos2 <- cbind(c(2, 5:7, 12), 1:5)
coefr3 <- matrix(0, 12, 7)
pos3  <- cbind(c(1:2, 5:7, 9, 12), 1:7)
coefr4 <- matrix(0, 12, 7)
pos4 <- cbind(c(1:2, 5, 6, 7, 9, 12), 1:7)
coefr5 <- matrix(0, 12, 5)
pos5 <- cbind(c(1, 2, 5, 7, 12), 1:5)
coefr6 <- matrix(0, 12, 5)
pos6 <- cbind(c(1, 2, 7, 9, 12), 1:5)
coefr1[pos1] <- 1
coefr2[pos2] <- 1
coefr3[pos3] <- 1
coefr4[pos4] <- 1
coefr5[pos5] <- 1
coefr6[pos6] <- 1

int <- 1:6
restrict.regMat <- c(list(coefr1, coefr2, coefr3, coefr4,
                          coefr5, coefr6, rep(list(NULL), 3)))

fit1 <- MultiChainLadder(da1, "SUR", int = int,
			restrict.regMat = restrict.regMat, model = "GMCL")
fit2 <- MultiChainLadder(da2, "OLS")
fit <- Join2Fits(fit1, fit2)
pred <- predict(fit)
mse <- Mse(fit, pred)
fit.sel <- JoinFitMse(fit, mse)

# 6. Munich Chain Ladder
fit.mucl <- MunichChainLadder(auto[[1]], auto[[2]])


### combine all the paid-to-incurred ratios

u1 <- summary(fit.scl)@Ultimate
r1 <- u1[, 1]/u1[, 2]
u2 <- summary(fit.mcl)@Ultimate
r2 <- u2[, 1]/u2[, 2]
u3 <- summary(fit.int)@Ultimate
r3 <- u3[, 1]/u3[, 2]
u4 <- summary(fit.full)@Ultimate
r4 <- u4[, 1]/u4[, 2]
u5 <- summary(fit.sel)@Ultimate
r5 <- u5[, 1]/u5[, 2]
r6 <- summary(fit.mucl)[[1]][, 6]
r6 <- c(r6, summary(fit.mucl)[[2]][2, 3])

ratios <- cbind(r1, r2, r3, r6, r4, r5)
dimnames(ratios)[[2]] <- c("SCL", "MCL", "GMCL1", "MuCL", "GMCL2", "GMCL3" )
ratios <- format(round(ratios*100 ,2), big.mark=",", scientific = FALSE) 
print(ratios, quote = FALSE)


### Estimated parameters
fit.sel@coefficients
fit.sel@residCov
summary(fit.sel)@residCor

## summary statistics
summary(fit.sel, portfolio = "1+3")@report.summary[[4]]


### residual plots
par(mfrow = c(2, 3))
plot(fit.scl, which.plot = 3:4)
plot(fit.mcl, which.plot = 3:4)
plot(fit.int, which.plot = 3:4)
plot(fit.full, which.plot = 3:4)
plot(fit.sel, which.plot = 3:4)
	
# histogram plots
r <- summary(fit.sel)@rstandard
par(mfrow = c(1,3))
for (i in 1:3){
	hist(r[, i], 20, freq = FALSE, xlim = c(-2, 2),
			 xlab = "Standardized Residuals",
			 main = paste("Histogram for Triangle", i))
	lines(density(r[, i]))
}