#
# Code to fit CNB model to a loss triangle by maximum likelihood
# Reference "Estimating Predictive Distributions for Loss Reserve Models"
#                            By Glenn Meyers
# Language - R (www.rproject.org)
#
# This code reads a loss development triangle from a comma seperated value file
# and (1) estimates loss development factors and the elr
#     (2) calculates a predictive mean and standard deviation for unpaid losses
#     (3) calculates a vector of probability densities and
#     (4) plots the densities
#
# Note - The cnb maximum likelihood search can take some time to converge.
#
# The data accessed by this program is illustrative only.
#
###################### Begin Functions #########################################
#
# Convert indata into a loss development rectangle for a given group
#
rectangle<-function(triangle){
    temp<-data.matrix(triangle)
    premium<-rep(0,55)
    ay<-rep(0,55)
    lag<-rep(0,55)
    loss<-rep(0,55)
    rrow<-0
    for (i in 1:9){
        rrow<-rrow+1
        ay[rrow]<-temp[i,1]-min(temp[,1])+1
        premium[rrow]<-temp[i,2]
        lag[rrow]<-1
        loss[rrow]<-temp[i,3]
        for (j in 2:(11-i)){
            rrow<-rrow+1
            premium[rrow]<-temp[i,2]
            ay[rrow]<-temp[i,1]-min(temp[,1])+1
            lag[rrow]<-j
            loss[rrow]<-temp[i,2+j]-temp[i,1+j]
            } #end j loop
        } #end i loop
    rrow<-rrow+1
    ay[rrow]<-temp[10,1]-min(temp[,1])+1
    premium[rrow]<-temp[10,2]
    lag[rrow]<-1
    loss[rrow]<-temp[10,3]
    rectframe<-data.frame(premium,ay,lag,loss)
    return(rectframe)
    } # end rectangle function
#
# estimate the step size for discrete loss distribution
#
estimate.h<-function(indata,fftn,limit.divisors){
    h<-sum(data.matrix(indata)[,2])/2^fftn
    h<-min(subset(limit.divisors,limit.divisors>h))
    return(h)
    }
#
# Pareto limited average severity
#
pareto.las<-function(x,theta,alpha){
  las<- ifelse(alpha==1,-theta*log(theta/(x+theta)),
               theta/(alpha-1)*(1-(theta/(x+theta))^(alpha-1)))
  return(las)
  }
#
# discretized Pareto severity distribution (different than that used in paper)
#
discrete.pareto<-function(limit,theta,alpha,h,fftn){
    dpar<-rep(0,2^fftn)
    m<-ceiling(limit/h)
    x<-h*0:m
    las<-rep(0,m+1)
    for (i in 2:(m+1)) {
        las[i]<-pareto.las(x[i],theta,alpha)
        } #end i loop
    dpar[1]<-1-las[2]/h
    dpar[2:m]<-(2*las[2:m]-las[1:(m-1)]-las[3:(m+1)])/h
    dpar[m+1]<-1-sum(dpar[1:m])
    return(dpar)
    } # end discrete.pareto function
#

#
# transformed dev to dev parameters
# transformed to enforce constraints
#
tdev.to.dev<-function(tdev){
    dev<-rep(0,11)
    dev[1]<-exp(-tdev[1]^2)/2
    dev[2]<-dev[1]*(1+exp(-tdev[2]^2))
    dev[3]<-min((1-sum(dev[1:2])),dev[2])*exp(-tdev[3]^2)
    dev[4]<-min((1-sum(dev[1:3])),dev[3])*exp(-tdev[4]^2)
    dev[5]<-min((1-sum(dev[1:4])),dev[4])*exp(-tdev[5]^2)
    dev[6]<-min((1-sum(dev[1:5])),dev[5])*exp(-tdev[6]^2)
    dev[7]<-min((1-sum(dev[1:6])),dev[6])*exp(-tdev[7]^2)
    dev[8]<-min((1-sum(dev[1:7])),dev[7])*exp(-tdev[8]^2)
    dev[9]<-min((1-sum(dev[1:8])),dev[8])*exp(-tdev[8]^2)
    dev[10]<-min((1-sum(dev[1:9])),dev[9])*exp(-tdev[8]^2)
    dev[11]<-tdev[9]^2
    temp<-sum(dev[1:10])
    dev[1:10]<-dev[1:10]/temp
    dev[11]<-temp*dev[11]
    return(dev)
    } # end tdev.to.dev function
#
# overdispersed Poisson loglikelihood (no sigma)  - from Dave Clark's paper
#
odpoisson.llike1<-function(par){
    dev<-tdev.to.dev(par)
    cdev<-c(dev[1:10],dev[1:9],dev[1:8],dev[1:7],dev[1:6],
          dev[1:5],dev[1:4],dev[1:3],dev[1:2],dev[1:1])
    elr<-dev[11]
    mu<-rdata$premium*cdev*elr
    loss<-rdata$loss
    llike1<-sum(loss*log(mu)-mu)
    return(llike1)
    } # end odpoisson.llike1
#
# probabililty of compound negative binomial distribution
#
dcnb.phiz<-function(x,eloss,c,fftn,h,esev,phiz){
    ecount<-eloss/esev
    phix<-(1-c*ecount*(phiz-1))^(-1/c)
    pv<-zapsmall(Re(fft(phix,inverse=TRUE))/2^fftn,digits=12)
    ind<-round(x/h)+1
    dcnb<-pv[ind]
    return(dcnb)
    } # end dcnb function
#
# compound negative binomial likelihood
#
cnb.llike1<-function(par){
    llike<-1:55
    dev<-tdev.to.dev(par)
    cdev<-c(dev[1:10],dev[1:9],dev[1:8],dev[1:7],dev[1:6],
          dev[1:5],dev[1:4],dev[1:3],dev[1:2],dev[1:1])
    elr<-dev[11]
    for (i in 1:55) {
        k<-rdata$lag[i]
        loss<-rdata$loss[i]
        eloss<-rdata$premium[i]*cdev[i]*elr
        llike[i]<-log(dcnb.phiz(loss,eloss,c,fftn,h,esev[k],phiz[,k]))
        }
    llike1<-sum(llike)
    return(llike1)
    } # end dcnb.llike1
#
# total reserve predictive distribution
#
totres.predictive.dcnb<-function(premium,dev,elr,esev,dp,c,fftn){
    phix<-complex(2^fftn,1,0)
    for (j in 2:10){
      dpp<-rep(0,2^fftn)
      lamp<-0
      for (k in (12-j):10) {
        eloss<-premium[j]*elr*dev[k]
        lam<-eloss/esev[k]
        lamp<-lamp+lam
        dpp<-dpp+dp[,k]*lam
        } # end k (dev) loop
      dpp<-dpp/lamp
      phix<-phix*(1-c*lamp*(fft(dpp)-1))^(-1/c)
      }  # end j (AY) loop
    pred<-zapsmall(Re(fft(phix,inverse=TRUE))/2^fftn,digits=12)
    return(pred)
    } #end totres predictive density
#
# total reserve mean and standard deviation
#
totres.predictive.mean.and.std<-function(premium,dev,elr,esev,dp,c,fftn,h){
    mean<-0
    var<-0
    for (j in 2:10){
      dpp<-rep(0,2^fftn)
      lamp<-0
      mp<-0
      for (k in (12-j):10) {
        eloss<-premium[j]*elr*dev[k]
        lam<-eloss/esev[k]
        lamp<-lamp+lam
        dpp<-dpp+dp[,k]*lam
        mp<-mp+eloss
        } #end k (dev) loop
        mean<-mean+mp
        var<-var+h*h*sum(dpp*0:(2^fftn-1)*0:(2^fftn-1))+c*mp^2
      } # end j (AY) loop
    std<-sqrt(var)
    return(c(mean,std,std/mean))
    } #end totres predictive mean and standard deviation function
#
# total reserve analysis
#
totres.Plot.Pred.Density<-function(premium,dev,elr,esev,dp,c,fftn,h,low=.5,hi=1.5){
    moments<-totres.predictive.mean.and.std(premium,dev,elr,esev,dp,c,fftn,h)
    rangemin<-low*moments[1]
    rangemax<-hi*moments[1]
    Range<-rangemin+(rangemax-rangemin)/200*0:200
    pred.density<-totres.predictive.dcnb(premium,dev,elr,esev,dp,c,fftn)
    Predictive.Density<-pred.density[(round(Range/h)+1)]
    par(ask=T)
    plot(Range/1000,Predictive.Density,xlab="Range (000,000)",type="l",yaxt="n")
    return(round(moments,3))
    } # End plot function
#
# End of functions
################################################################################
# Begin main program
#
# read in data from a csv file
#
indata<-read.csv("Demo Triangle.csv")
#
# set up discretized claim severity distribution
#
theta<-c(10,25,50,75,100,125,150,150,150,150)
alpha<-rep(2,10)
limit<-1000
limit.divisors<-c(5,10,20,25,40,50,100,125,200,250,500,1000)
fftn<-14
h<-estimate.h(indata,fftn,limit.divisors)
phiz<-matrix(0,2^fftn,10)
dp<-matrix(0,2^fftn,10)
esev<-rep(0,10)
for (k in 1:10){
    dp[,k]<-discrete.pareto(limit,theta[k],alpha[k],h,fftn)
    esev[k]<-pareto.las(limit,theta[k],alpha[k])
    phiz[,k]<-fft(dp[,k])
    } #end k loop
c<-.01

############ execute these lines to use the Clark log likelihood #######
par<-rep(1,9)
rdata<-rectangle(indata)
ml.odp<-optim(par,odpoisson.llike1,control=list(fnscale=-1,maxit=10000))
dev<-tdev.to.dev(ml.odp$par)
dev.odp<-dev[1:10]
elr.odp<-dev[11]
############ execute these lines to use the cnb log likelihood #######
init.time<-proc.time()[3]/60
ml.cnb<-optim(ml.odp$par,cnb.llike1,control=list(fnscale=-1,maxit=2000))
dev<-tdev.to.dev(ml.cnb$par)
dev.cnb<-dev[1:10]
elr.cnb<-dev[11]
############  print sample results ###################################
run.time<-proc.time()[3]/60-init.time
run.time
dev.odp
dev.cnb
dev.odp/dev.cnb
elr.odp
elr.cnb
elr.odp/elr.cnb
premium<-data.matrix(indata)[,2]
totres.Plot.Pred.Density(premium,dev.cnb,elr.cnb,esev,dp,c,fftn,h)
totres.predictive.mean.and.std(premium,dev.cnb,elr.cnb,esev,dp,c,fftn,h)
############ calculate mean and standard deviation from the density function
temp<-totres.predictive.dcnb(premium,dev.cnb,elr.cnb,esev,dp,c,fftn)
m1<-h*sum(0:(2^fftn-1)*temp)
s1<-sqrt(h*h*sum(0:(2^fftn-1)*0:(2^fftn-1)*temp)-m1^2)
m1
s1