Reserving_based_on_log-incremental_payments.R

## Copyright Markus Gesmann, January 2013
## See: http://lamages.blogspot.co.uk/2013/01/reserving-based-on-log-incremental.html
## Also: www.actuaries.org.uk/system/files/documents/pdf/crm2-D5.pdf
## 
## Incremental claims triangle
tri <- t(matrix(
  c(11073,  6427, 1839, 766,
    14799,  9357, 2344,  NA,
    15636, 10523,   NA,  NA,
    16913,    NA,   NA,  NA), 
  nc=4, dimnames=list(origin=0:3, dev=0:3)))

m <- dim(tri)[1]; n <- dim(tri)[2]
dat <- data.frame(
  origin=rep(0:(m-1), n),
  dev=rep(0:(n-1), each=m),
  value=as.vector(tri))

## Add dimensions as factors
dat <- with(dat, data.frame(origin, dev, cal=origin+dev, 
                            value, logvalue=log(value),
                            originf=factor(origin),
                            devf=as.factor(dev),
                            calf=as.factor(origin+dev)))

rownames(dat) <- with(dat, paste(origin, dev, sep="-"))
dat <- dat[order(dat$origin),]
dat


Fit <- lm(logvalue ~ originf + devf + 0, data=dat)
summary(Fit)

# Resdiual plots
op <- par(mfrow=c(2,2))
attach(model.frame(Fit))
plot.default(rstandard(Fit) ~ originf,
             main="Residuals vs. origin years")
abline(h=0, lty=2)
plot.default(rstandard(Fit) ~ devf,
             main="Residuals vs. dev. years")
abline(h=0, lty=2)
with(na.omit(dat), 
     plot.default(rstandard(Fit) ~ calf,
                  main="Residuals vs. payments years"))
abline(h=0, lty=2)
plot.default(rstandard(Fit) ~ logvalue,
             main="Residuals vs. fitted")
abline(h=0, lty=2)
detach(model.frame(Fit))
par(op)

# Model design matrix
dm <- model.matrix(formula(Fit), dat=model.frame(Fit))
dm

## Future design matrix
fdm <- model.matrix( ~ originf + devf + 0, data=dat[is.na(dat$value),])
fdm

varcovar <- fdm %*% vcov(Fit) %*% t(fdm)
## fdm %*% solve(t(dm)%*%dm) %*% t(fdm) *sigma^2, or shorter
round(varcovar, 4)

sigma <- summary(Fit)$sigma
sigma

Var <- varcovar + sigma^2
VarY <- Var[row(Var)==col(Var)]
Y <- fdm %*% coef(Fit) 
# predict(Fit, newdata=dat[is.na(dat$value), c("originf", "devf")])
P <- exp(Y + VarY/2)
VarP <- exp(2*Y + VarY)*(exp(VarY)-1)
seP <- sqrt(VarP)
i <- fdm %*% c((1:m)-1, rep(0, (n-1)))
j <- fdm %*% c(rep(0, (m-1)), (1:n)-1)         
Results <- data.frame(i,j, Y, VarY, P, VarP, seP)
Results ## Page D5.13

newData <- dat[is.na(dat$logvalue), c("originf", "devf")]
Pred <- predict(Fit, newdata=newData, se.fit=TRUE)
Y <- Pred$fit
VarY <- Pred$se.fit^2 + Pred$residual.scale^2 
fdm <- model.matrix(~ originf + devf + 0, data=newData)

lower.tri <- xtabs(P ~ i+j, dat=Results)

Full.Incr.Triangle <- tri
Full.Incr.Triangle[row(tri) > (nrow(tri) + 1 - col(tri))] <- 
  lower.tri[row(lower.tri) > (nrow(lower.tri) - col(lower.tri))]

Full.Cum.Triangle <- apply(Full.Incr.Triangle, 1, cumsum) 
Full.Cum.Triangle

CoVar <-  sweep(sweep((exp(varcovar)-1) , 1, P, "*"), 2, P, "*")
CoVar[col(CoVar)==row(CoVar)] <- 0
round(CoVar,0) 
OverallVar <- sum(CoVar) + sum(VarP)
Total.SE <- sqrt(OverallVar)
Total.SE

Overall.Reserve <- sum(lower.tri)
Overall.Reserve

Total.SE / Overall.Reserve

## Compare to Mack-Chain-Ladder
library(ChainLadder)
M <- MackChainLadder(incr2cum(tri), est.sigma="Mack")
M


## Second example

dat <- data.frame(
  origin=rep(0:6, each=7),  
  dev=rep(0:6, 7),
  value= c(3511,  3215,  2266,  1712,  1059,	587,	340,
           4001,	3702,	2278,	1180,	956,	629,	NA,
           4355,	3932,	1946,	1522,	1238,	NA,	NA,
           4295,	3455,	2023,	1320,	NA,	NA,	NA,
           4150,	3747,	2320,	NA,	NA,	NA,	NA,
           5102,	4548,	NA,	NA,	NA,	NA,	NA,
           6283,	NA,	NA,	NA,	NA,	NA,	NA))

dat <- with(dat, 
            data.frame(origin, dev, cal=origin+dev,
                       value, logvalue=log(value), 
                       a6 = ifelse(origin == 6, 1, 0),
                       a5 = ifelse(origin == 5, 1, 0),
                       d = ifelse(dev < 1, 1, 0),
                       s = ifelse(dev < 1, 0, dat$dev)))        

summary(Fit <- lm(logvalue ~ a5 + a6 + d + s, data=na.omit(dat)))