setwd("C:/working directory")
insurer.data="C:/CAS Loss Reserve Database/comauto_pos.csv"
grpcode="353"
losstype="incloss"  #"incloss" if incurred loss or "cpdloss" if paid loss
a=read.csv(insurer.data)
#
# function to get Schedule P triangle data given ins group and line of business
#
ins.line.data=function(g.code){
  b=subset(a,a$GRCODE==g.code)
  name=b$GRNAME
  grpcode=b$GRCODE
  w=b$AccidentYear
  d=b$DevelopmentLag
  cum_incloss=b[,6]
  cum_pdloss=b[,7]
  bulk_loss=b[,8]
  dir_premium=b[,9]
  ced_premium=b[,10]
  net_premium=b[,11]
  single=b[,12]
  posted_reserve97=b[,13]
  # get incremental paid losses - assume data is sorted by ay and lag
  inc_pdloss=numeric(0)
  for (i in unique(w)){
    s=(w==i)
    pl=c(0,cum_pdloss[s])
    ndev=length(pl)-1
    il=rep(0,ndev)
    for (j in 1:ndev){            
      il[j]=pl[j+1]-pl[j]
      }
    inc_pdloss=c(inc_pdloss,il)
    }
  data.out=data.frame(grpcode,w,d,net_premium,dir_premium,ced_premium,
     cum_pdloss,cum_incloss,bulk_loss,inc_pdloss,single,posted_reserve97)
  return(data.out)
  }
#
# read and aggregate the insurer data and 
# set up training and test data frames
cdata=ins.line.data(grpcode)
w=cdata$w-1987
d=cdata$d
o1=100*d+w
o=order(o1)
w=w[o]
d=d[o]
premium=cdata$net_premium[o]
cpdloss=cdata$cum_pdloss[o]
incloss=cdata$cum_incloss[o]-cdata$bulk_loss[o]
adata=data.frame(grpcode,w,d,premium,cpdloss,incloss)
rdata=subset(adata,(adata$w+adata$d)<12)
maxlevel=log(2*max(rdata$premium))
numw=length(unique(rdata$w))
if(losstype=="incloss") rloss=rdata$incloss else rloss=rdata$cpdloss
if(losstype=="incloss") aloss=adata$incloss else aloss=adata$cpdloss
#
# run Mack model using ChainLadder
#
library(ChainLadder)
rtriangle=as.triangle(rdata,origin="w",dev="d",value=losstype)     #######################
mcl=MackChainLadder(rtriangle,est.sigma="Mack")
print(rtriangle)
print(" ")
print(" ")
#
# run JAGS model
#
library(rjags)
jags=jags.model('LCL2-JAGS.txt',
                   data = list('logloss'= log(rloss),
                               'maxlev' = maxlevel,
                               'numlev' = numw,
                               'w'      = rdata$w,
                               'd'      = rdata$d),
                   inits = list('.RNG.name'= "base::Marsaglia-Multicarry",
                                '.RNG.seed'= 12345),
                   n.chains = 1,
                   n.adapt = 100000,
                   )
#
# get the output from JAGS
#
#if necessary do #
#update(jags,100000)  
jagout=coda.samples(jags,c("level","dev","sigd2","z"),100000,thin=10)
plot(jagout)
print(summary(jagout))
b=as.matrix(jagout)
dev=b[,1:10]
level=b[,11:20]
sigd2=b[,21:30]
z=b[,31]
#
# simulate loss statistics by accident year reflecting both parameter and total risk
# for the JAGS model
#
set.seed(12345)
sim.ayd10=matrix(0,dim(b)[1],10)
ss.ayd10=rep(0,10)
ms.ayd10=rep(0,10)  
w=1
for(i in 1:dim(b)[1]){
  sim.ayd10[i,w]=rlnorm(1,meanlog=level[i,w]+dev[i,10],sdlog=sqrt(sigd2[i,10]))
}
ss.ayd10[w]=sd(sim.ayd10[,w])
ms.ayd10[w]=mean(sim.ayd10[,w])
for (w in 2:10){
  for(i in 1:dim(b)[1]){
    sim.ayd10[i,w]=rlnorm(1,meanlog=level[i,w]+dev[i,10]+
      z[i]*(log(sim.ayd10[i,w-1])-level[i,w-1]-dev[i,10]),sdlog=sqrt(sigd2[i,10]))
  }
  ss.ayd10[w]=sd(sim.ayd10[,w])
  ms.ayd10[w]=mean(sim.ayd10[,w])
}
ms.td10=mean(rowSums(sim.ayd10[,2:10]))
ss.td10=sd(rowSums(sim.ayd10[,2:10]))
LCL.Estimate=round(ms.ayd10)
LCL.S.E=round(ss.ayd10)
LCL.CV=round(LCL.S.E/LCL.Estimate,4)
act=sum(subset(aloss,adata$d==10)[2:10])
pred.LCL=rowSums(sim.ayd10[,2:10])
pct.LCL=sum(pred.LCL<=act)/length(pred.LCL)*100
#
# get loss statistics by accident year reflecting total risk from Mack model
# from Mack model
#
Mack.Estimate=round(summary(mcl)$ByOrigin[,3])
Mack.S.E=round(summary(mcl)$ByOrigin[,5])
Mack.CV=round(Mack.S.E/Mack.Estimate,4)
Mack.total.ult=round(sum(summary(mcl)$ByOrigin[2:10,3]))
Mack.total.se=round(summary(mcl)$Totals[5,1])
Mack.total.cv=round(Mack.total.se/Mack.total.ult,4)
Mack.sig=sqrt(log(1+Mack.total.cv^2))
Mack.mu=log(Mack.total.ult)-Mack.sig^2/2
pct.Mack=round(plnorm(act,Mack.mu,Mack.sig)*100,2)

#
# put JAGS and Mack accident year statistics into a data frame
#
W=c(1:10,"Total w=2:10")
LCL.Estimate=c(LCL.Estimate,round(ms.td10))
LCL.S.E=c(LCL.S.E,round(ss.td10))
LCL.CV=c(LCL.CV,round(ss.td10/ms.td10,4))
Mack.Estimate=c(Mack.Estimate,Mack.total.ult)
Mack.S.E=c(Mack.S.E,Mack.total.se)
Mack.CV=c(Mack.CV,Mack.total.cv)
risk=data.frame(W,LCL.Estimate,LCL.S.E,LCL.CV,Mack.Estimate,
                Mack.S.E,Mack.CV)
#
# create a histogram of the correlation coefficient z
#
par(mfrow=c(2,2))
hist(z, main="Histogram of the Coefficient of Correlation",xlab="")
par(mfrow=c(1,1))
#
# graphic retrospective test for ten year cumulative loss by accident year in JAGS model
#
par(mfrow=c(1,1))
samp=sample(1:dim(b)[1],100)
act10d=subset(aloss,adata$d==10)
yylim=range(act10d,sim.ayd10[samp,10])
LCL.w=1:10+0.1
Mack.w=1:10-0.1
xxlim=range(LCL.w,Mack.w)
Mack.sim=matrix(act10d[1],length(samp),10)
Mack.sim.sig=sqrt(log(1+Mack.CV[1:10]^2))
Mack.sim.mu=log(Mack.Estimate[1:10])-Mack.sim.sig^2/2
for (w in 2:10){
  for(i in 1:length(samp)){
    Mack.sim[i,w]=rlnorm(1,meanlog=Mack.sim.mu[w],sdlog=Mack.sim.sig[w])
  }
}
yylim=range(act10d,sim.ayd10[samp,10],Mack.sim)
plot(1:10,act10d,main="Retrospective Test of Loss Intervals",type="n",ylim=yylim,xlim=xxlim,  
     xlab="Accident Year",ylab="Cumulative Loss at 10 Years")
par(new=T)
plot(1:10,act10d,type="b",lwd=3,col="blue",main="",xlab="",ylab="",ylim=yylim,xlim=xxlim)
for(i in 1:length(samp)){
  par(new=T)
  plot(LCL.w,sim.ayd10[samp[i],],main="",xlab="",ylab="",type="p",pch=20,ylim=yylim,xlim=xxlim)
  par(new=T)
  plot(Mack.w,Mack.sim[i,],main="",xlab="",ylab="",type="p",pch=20,ylim=yylim,col="gray",xlim=xxlim)
}
#
# print results
#
Actual=c(act10d,act)
risk=cbind(risk,Actual)
print(risk)
write.csv(risk,file="risk 2.csv",row.names=F)
print(paste("LCL Percentile of Total Loss=",pct.LCL))
print(paste("Mack Percentile =",pct.Mack))