##===========================================================
 ## Chain Ladder Linear Regression
 ## Markus Gesmannn, July 2006, 4. Septemper 2006
 ## Collection of ideas on how to use linear regression for reserving,
 ## some were taken from the following paper:
 ## http://www.casact.org/pubs/proceed/proceed00/00245.pdf
 ## In order to run the examples below, you need my reserving package
 ## available from the toolkit wiki:
 ## http://toolkit.pbwiki.com/f/reserving_0.1-2.zip 
 ## The idea behind this concept is to use the build in functionality of R
 ## to check and modify and plot models and residuals
 ## Copy and paste this code into R
 ##===========================================================


 lmChainLadder <- function(data, formula='y~x+0',
                          weights.FUN=function(x){1/abs(x)},
                          ...){
 ## Think of x= claims at dev. period i, and 
 ## y=claims at dev. period i+1
 ## The ordinary Chain Ladder ratio is nothing but 
 ## a weighted linear regression through the origin.	
 ## The data should be named "YOA", "Dev", "Amount"
 ## and in a table structure, not a triangle
 ##
  cl <- match.call()
  mf <- match.call(expand.dots = FALSE)
  m <- match(c("data", "formula", "weights.FUN"), names(mf), 0)
  mf <- mf[c(1, m)]
 ##
 ##
  data <- na.omit(data)
 ## 
  spL=split(data, data$Dev) 
  dev.period=names(spL)
  chain.ladder.model <- lapply(dev.period, function(x){ 
    i=as.numeric(as.character(x))
    if(i<length(dev.period)){     
      y <- spL[[as.character(i+1)]]$Amount
      yn <- length(y)
      x <- spL[[as.character(i)]]$Amount[1:yn]
      if(length(x)>1){
        model <- lm(as.formula(formula), data=data.frame(x=x, y=y), 
        	weights=weights.FUN(x),...)
      }else{
        model <- lm(y~x+0, data=data.frame(x=x, y=y), weights=weights.FUN(x),...)
      }
      return(model)
    }
   }
  )
  names(chain.ladder.model) <- dev.period
  Result <- list()
  Result[["Call"]] <- mf
  Result[["model"]] <- chain.ladder.model
  Result[["formula"]] <- as.formula(formula)
  Result[["weights.FUN"]] <- weights.FUN
  Result[["data"]] <- data
  class(Result) <- c("lmChainLadder", class(Result))
  return(Result)
 ##
 }

 ##===========================================================

 residuals.lmChainLadder <- function(object){
 ## function to extract the residuals	
  model <- object$model
  data <- object$data
  f <- names(model)
  f <- f[-length(f)]
 ## 
  .Res <- do.call("c",
                  lapply(f, function(x){
                    .wghts <- weights(model[[x]])
                    .wghts <- ifelse(is.null(.wghts),1,.wghts)
                    residuals(model[[x]])*sqrt(.wghts)/
                      summary.lm(model[[x]])$sigma
                  })
                  )
  ##
 .Fitted <- do.call("c",
                    lapply(f, function(x){
                      fitted(model[[x]])
                    })
                    )
  .YOA <- data[,1]
  if(class(.YOA)=="factor"){
    .YOA <- as.numeric(levels(.YOA))[.YOA]
  }
  .YOA <- .YOA[.YOA!=min(.YOA)]
  .Dev <- data[,2]
  if(class(.Dev)=="factor"){
    .Dev <- as.numeric(levels(.Dev))[.Dev]
  }
  .Dev <- .Dev[.Dev!=min(.Dev)]
  ##
  Result <- data.frame(YOA=.YOA, Dev=.Dev, Calendar=.YOA+.Dev-1,
                       Residuals=.Res, Fitted=.Fitted)
  class(Result) <- c("ChainLadderResiduals", class(Result))
  return(Result)
 }

 ##===========================================================

 predict.lmChainLadder <- function(object){
 ## function to predict the ultimates	
  data <- object$data
  model <- object$model
  pred.triangle <- Long.Triangle(data)[,-1]
  if(!match(names(pred.triangle)[1],c("UWY", "YOA", "YoA"),FALSE)){
    pred.tirangle <- pred.triangle[,-1]
  }
  for(i in as.numeric(names(model))){
    if(i<length(names(model))){
      fp <- which( is.na(pred.triangle[,(i+1)]) & !is.na(pred.triangle[,i])  )
      if(length(fp)>0){
        pred.triangle[fp,(i+1)] <- predict(model[[as.character(i)]],
                                           newdata=data.frame(x=pred.triangle[fp,i]))
      }
    }
  }
  return(pred.triangle)    
 }

 ##===========================================================

 plot.ChainLadderResiduals <- function(x,
                                      FUN=function(x){x/sd(x)},
                                      na.rm=TRUE
                                      ){
  if(na.rm){
    x <- na.omit(x)
  }
  .op <- par(mfrow=c(2,2))
  plot(FUN(Residuals) ~ Dev, data=x)
  abline(h=0)
  lines(lowess(y=(FUN(x$Residuals)), x=x$Dev), col=2)
  ## 
  plot(FUN(Residuals) ~ YOA, data=x)
  abline(h=0)
  lines(lowess(y=(FUN(x$Residuals)),x=x$YOA), col=2)
  ##
  plot(FUN(Residuals) ~ Calendar, data=x)
  abline(h=0)
  lines(lowess(y=(FUN(x$Residuals)), x=x$Calendar), col=2)
  ##
  plot(FUN(Residuals) ~ Fitted, data=x)
  abline(h=0)
  lines(lowess(y=(FUN(x$Residuals)), x=x$Fitted), col=2)
  ##
  par(oma=c(0,0,2,0))
   # question to R-help if the function extraction is correct"
  .funbody <- as.character(body(FUN))
  .funtext <- ""
  for(i in seq(along=.funbody)){
    if(.funbody[i]!="{")
      .funtext <- paste(.funtext,.funbody[i])
  }
  title(paste("Residual Plot, FUN=",.funtext,sep="") , outer=TRUE) 
  par(.op)
 }

 ##===========================================================

 print.lmChainLadder <- function(x,...){
  cat("Call: ")
  print(x$Call)
  cat("n")
  model <- x$model
  f <- names(model)
  f <- f[-length(f)]
  result <- do.call("rbind",
                 lapply(f, function(x){                  
                   cbind(summary.lm(model[[x]])$coef,
                         summary.lm((model[[x]]))$sigma,
                         summary.lm((model[[x]]))$r.squared) 
                 })
                 )
  colnames(result) <- c("Estimate Std.", "Error", "t value", "Pr(>|t|", "Sigma", "R-squared")
  rownames(result) <- paste(rownames(result))
  print(result)
 }

 ##===========================================================

 coef.lmChainLadder <- function(object){ 
  model <- object$model
  f <- names(model)
  f <- f[-length(f)]
  result <- lapply(f, function(x){                  
                   c(coef(model[[x]]) )
                 })
  names(result) <- f
  return(result)
 }

 ##===========================================================
 ## Examples

  require(reserving)
  data(GenIns)
  par(ask=TRUE)
  print("Data:")
  print(GenIns)
  # convert triangle into table
  L <- Triangle.Long(GenIns)
  cat("Ordinary Chain Laddern")
  M <- lmChainLadder(L,formula='y~x+0', weights.FUN=function(x){1/abs(x)})
  print(M)
  plot(residuals(M))
  pM <- predict(M)
  print(pM) 
  cat("Regression without weighting and the restriction to the interceptn")
  G <- lmChainLadder(L,formula='y~x', weights.FUN=function(x){rep(1,length(x))})
  print(G)
  plot(residuals(G))
  pG <- predict(G)
  print(pG)
  ##
  cat("compare the prediction of the above modelsn")
  print(pM/pG)
  matplot(t(pM/pG), t="o")
  ##
  cat("Zenwirth's approachn")
  IncGenIns <- incTria(GenIns)
  ## beawre when you have negative incrementals
  LogIncGenIns <- log(IncGenIns)
  LogIncL <- Triangle.Long(LogIncGenIns)
  Z <- lmChainLadder(LogIncL,
                     formula='y~x',
                     weights.FUN=function(x){rep(1,length(x))})
  ##
  print(Z)
  plot(residuals(Z))
  pZ <- cumTria(exp((predict(Z))))
  print(pZ)
  ##
  cat("Compare Zehnwirth's approach with ordiniary Chain Laddern")
  print(pZ/pM)
  matplot(t(pM/pZ), t="o")
  par(ask=FALSE)