#
# 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)
  }
#
# estimate the step size for discrete loss distribution
#
estimate.h<-function(prob.max.loss,fftn){
    h<-prob.max.loss/2^fftn
    limit.divisors<-c(.5,1,2.5,5,10,20,25,40,50,100,125,200,250,500,1000,2500)
    h<-min(subset(limit.divisors,limit.divisors>h))
    return(h)
    }
#
# discretized Pareto severity distribution
#
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
#
# Begin main program #########################################################
#
# Input parameters
#
theta<-10
alpha<-2
limit<-1000
fftn<-14
c<-.01
pml.multiplier=10
expected.loss=10000
#
# Set up discretized severity distributions
#
h=estimate.h(pml.multiplier*expected.loss,fftn)
dp<-discrete.pareto(limit,theta,alpha,h,fftn)
#
# Calculate probabilities for aggregate loss distribution
#
esev<-pareto.las(limit,theta,alpha)
lambda=expected.loss/esev
phix<-(1-c*lambda*(fft(dp)-1))^(-1/c)
agg.prob=zapsmall(Re(fft(phix,inverse=TRUE))/2^fftn,digits=12)
#
# Check for fft recycling - check.expected.loss should be equal to expected loss
#
x=h*0:(2^fftn-1)
check.expected.loss=sum(x*agg.prob)
check.expected.loss
#
# Plot density of aggregate loss distribution
#
windows(record=T)
plot(x,agg.prob)
cum.prob=cumsum(agg.prob)
trim.tails=cum.prob<.999999
trim.tails=(cum.prob>.000001)&trim.tails
plot(x[trim.tails],agg.prob[trim.tails],main="Collective Risk Model",
     type="l",lwd=3,col="red")
#
# Collective risk model simulation
#
nsim=10000
agg.loss=rep(0,nsim)
set.seed(1234)
for (i in 1:nsim){
  chi=rgamma(1,1/c,1/c) # random gamma variate with E[chi]=1 and Var[chi]=c
  n=rpois(1,chi*lambda)
  p=runif(n)
  agg.loss[i]=sum(pmin(theta*(1/(1-p)^(1/alpha)-1),limit))
  }
#
# Plot histogram of agg.loss
#
par(new=T) # Put new histogram on top of density plot
hist(agg.loss,xlim=range(x[trim.tails]),axes=F,main="",ylab="",xlab="")
legend("topright",c("FFT Inversion","Simulation"),col=c("red","black"),
        lwd=c(3,1))