Reserving_based_on_log-incremental_payments2.R

# Markus Gesmann, January 2013
# http://lamages.blogspot.co.uk/2013/01/reserving-based-on-log-incremental_15.html
# Page D5.17
tri <- t(matrix(
  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), nc=7))

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)))

op <- par(mfrow=c(2,1), mar=c(4,4,2,2))
with(dat, interaction.plot(x.factor=dev, trace.factor=origin, 
                           response=value))
points(dat$devf, dat$value, pch=16, cex=0.5)
with(dat, interaction.plot(x.factor=dev, trace.factor=origin, 
                           response=logvalue))
points(dat$devf, dat$logvalue, pch=16, cex=0.5)
par(op)

dat <- with(dat, 
            data.frame(dat,
                       a6 = ifelse(origin == 6, 1, 0),
                       a5 = ifelse(origin == 5, 1, 0),
                       d = ifelse(dev < 1, 1, 0),
                       s = ifelse(dev < 1, 0, dev)))        

dat <- dat[order(dat$origin),] 

# Page D5.20/21
summary(Fit1 <- lm(logvalue ~ 0 + originf + d + s, data=na.omit(dat)))

# Page D5.21/22
summary(Fit2 <- lm(logvalue ~ originf + d + s, data=na.omit(dat)))

# Page D5.23
summary(Fit3 <- lm(logvalue ~ a5 + a6 + d + s, data=na.omit(dat)))

# Resdiual plots
op <- par(mfrow=c(2,2), mar=c(4,4,2,2))
attach(model.frame(Fit3))
with(na.omit(dat), 
     plot.default(rstandard(Fit3) ~ origin,
                  main="Residuals vs. origin years"))
abline(h=0, lty=2)

with(na.omit(dat), 
     plot.default(rstandard(Fit3) ~ dev,
                  main="Residuals vs. dev. years"))
abline(h=0, lty=2)
with(na.omit(dat), 
     plot.default(rstandard(Fit3) ~ cal,
                  main="Residuals vs. payments years"))
abline(h=0, lty=2)
plot.default(rstandard(Fit3) ~ logvalue,
             main="Residuals vs. fitted")
abline(h=0, lty=2)
detach(model.frame(Fit3))
par(op)

log.incr.predict <- function(
  model, # lm output
  newdata, # same argument as in predict
  origin.var="origin", # name of the origin column in newdata
  dev.var="dev" # name of the dev column in newdata
){
  origin <- newdata[[origin.var]] 
  dev <- newdata[[dev.var]]
  
  Pred <- predict(model, newdata=newdata, se.fit=TRUE)
  
  Y <- Pred$fit
  VarY <- Pred$se.fit^2 + Pred$residual.scale^2 
  P <- exp(Y + VarY/2)
  VarP <-  exp(2*Y + VarY)*(exp(VarY)-1)
  seP <- sqrt(VarP)
  
  ## Recreate formula to derive the model.frame and future design matrix   
  model.formula <- as.formula(paste("~", as.character(formula(model)[3])))
  ## See also package formula.tool
  
  mframe <- model.frame(model.formula, data=newdata)
  fdm <- model.matrix(model.formula, data=newdata)
  varcovar <- fdm %*% vcov(model) %*% t(fdm)
  OverallVar <- sqrt( t(P) %*% (exp(varcovar)-1) %*% P )
  Total.SE <- sqrt(OverallVar)
  Total.Reserve <- sum(P)
  
  # Prepare output
  Incr=data.frame(origin, dev, Y, VarY, P, seP, CV=seP/P)
  out <- list(Forecast=Incr[order(newdata[[origin.var]]),],              
              Totals=data.frame(Total.Reserve, 
                                Total.SE=Total.SE, 
                                CV=Total.SE/Total.Reserve)
  )
  
  return(out)
}

tail.years <- 6
# Create a data frame for the future periods
fdat <- data.frame(
  origin=rep(0:(m-1), n+tail.years),
  dev=rep(0:(n+tail.years-1), each=m)
)
fdat$cal <- with(fdat, origin + dev) 

fdat <- with(fdat, 
             data.frame(fdat,
                        originf = factor(origin),
                        a6 = ifelse(origin == 6, 1, 0),
                        a5 = ifelse(origin == 5, 1, 0),
                        d = ifelse(dev < 1, 1, 0),
                        s = ifelse(dev < 1, 0, dev)))  

FM2 <- log.incr.predict(Fit2, subset(fdat, cal>6))
FM2 # Page D5.32/33

FM3 <- log.incr.predict(Fit3, subset(fdat, cal>6))
FM3 # Page D5.35

library(ChainLadder)
M <- MackChainLadder(incr2cum(tri), est.sigma="Mack", 
                     tail=1.05, tail.se=0.02)
M