#
# overdispersed Poisson loglikelihood - from Dave Clark's paper
#
beta.odpoisson.llike<-function(par){
    llike=1:55
    cdev1=pbeta((1:10)/10,par[1],par[2])
    dev=cdev1-c(0,cdev1[1:9])
    idev<-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<-par[3]
    mu<-premium*idev*elr
    llike1<-sum(loss*log(mu)-mu)
    return(llike1)
    } # end beta.odpoisson.llike
#
# Begin main program
#
# Set up initial parameters
#
premium=rep(50000,55)
ay=c(rep(1,10),rep(2,9),rep(3,8),rep(4,7),rep(5,6),rep(6,5),
     rep(7,4),8,8,8,9,9,10)
lag=c(1:10,1:9,1:8,1:7,1:6,1:5,1:4,1:3,1,2,1)
base.a=1.5
base.b=5
base.elr=.7
#
# Calculate expected development factors and losses
#
cum.dev=pbeta((1:10)/10,base.a,base.b)
inc.dev=cum.dev-c(0,cum.dev[1:9])
e.loss=1:55
for (i in 1:55) {
  e.loss[i]=base.elr*premium[i]*inc.dev[lag[i]]
  }
windows(record=T)
plot(1:10,inc.dev,type="l",lwd=3,col="red",
     main="Incremental Development Factors")
#
# Simulate losses from the overdispersed Poisson distribution and
# (1) estimate the parameters of the model by maximum likelihood
# (2) plot the parameters to see parameter risk
#
nsim=100
a=1:nsim
b=1:nsim
elr=1:nsim
set.seed(1234)
sigma.sq=.1*max(e.loss) ######### This deserves scrutiny! ################
lambda=e.loss/sigma.sq
for (i in 1:nsim) {
  x=rpois(55,lambda)
  loss=x*sigma.sq
  ml<-optim(c(base.a,base.b,base.elr),beta.odpoisson.llike,
            control=list(fnscale=-1,maxit=10000))
  a[i]=ml$par[1]
  b[i]=ml$par[2]
  elr[i]=ml$par[3]
  cdev=pbeta((1:10)/10,a[i],b[i])
  idev=cdev-c(0,cdev[1:9])
  par(new=T)
  plot(1:10,idev,type="l",axes=F,xlab="",ylab="",ylim=range(inc.dev))
  }
#
par(new=T)
plot(1:10,inc.dev,type="l",lwd=3,col="red",ylim=range(inc.dev),axes=F,main="")
legend("topright",legend=c("True","Simulated"),lwd=c(3,1),col=c("red","black"))
#
plot(a,b)
hist(elr)