library(MASS)
library(stats4)
#### Get the parameters for frequency and severity distributions
n <- 1000000 ### 1 million scenarios
lambda <- 3 ### Poisson rate
lognor.mu <- 7 ### Mu for lognormal distribution
lognor.sigma <- 1.3 ### Sigma for lognormal distribution
#### Misc. functions
### Function to read a vector of length 1+n and produce an n-dimensional vector
### with k 1's and n-k 0's
### where k = first element of the input vector (x[1])
Set.zeroes = function(x) {
multiplier <- c(rep(1,x[1]), rep(0,length(x) -1-x[1]))
x <- x[-1]*multiplier
}
### Monte Carlo Simulation for Gross Aggregate loss data
### Simulate n different claim counts from Poisson distribution assuming frequency ~ Poisson (lambda)
claim.count <- rpois(n, lambda)
max.claim.count <- max(claim.count) ### Builds a finite matrix of losses
### Create a matrix of n x max.claim.count individual losses
ind.losses <- rlnorm(n * max.claim.count, meanlog = lognor.mu, sdlog = lognor.sigma)
loss.amounts <- matrix(ind.losses, ncol = max.claim.count)
loss.matrix <- cbind(claim.count,loss.amounts)
head(loss.matrix)
## claim.count
## [1,] 7 2215.1474 366.1524 845.0304 1105.1524 963.3906
## [2,] 2 1021.0718 4350.5106 451.3270 534.3410 595.7512
## [3,] 4 504.5278 2259.4447 344.6350 1310.5767 399.2020
## [4,] 2 8230.7819 327.0207 7166.7646 5322.3117 10330.6050
## [5,] 1 682.2760 451.6263 724.5964 270.8565 636.3139
## [6,] 4 160.8754 531.7196 2644.2613 2808.4345 4834.0660
##
## [1,] 466.1302 102.3989 2918.1715 3533.0806 1227.48992 425.3582
## [2,] 1633.5624 148.2146 2529.0733 2474.5248 7962.19201 1659.8151
## [3,] 380.9781 2586.4789 595.8826 1866.1372 83.73633 1764.8210
## [4,] 4870.6072 8433.1841 2253.9975 998.0897 1349.67703 521.4459
## [5,] 3642.0152 1169.2948 13936.2683 1750.5243 10062.99699 124.1837
## [6,] 3296.2546 1254.2181 1233.0513 679.6782 779.43142 1585.0856
##
## [1,] 367.7808 5892.8245 17058.98034
## [2,] 262.9563 949.5997 453.50713
## [3,] 146.2620 1914.2644 97.45901
## [4,] 907.8022 1800.4514 640.99418
## [5,] 6915.4853 2030.2394 3832.15116
## [6,] 5900.3483 245.8644 620.99224
### For each scenario retain only the first claim.count losses
losses <- t(apply(loss.matrix, 1, Set.zeroes)) ### transpose the matrix
agg.losses <- apply(losses, 1, sum) ### calculate the aggregate losses for each scenario
### Visualize the agg.losses
library(ggplot2)
breaks <- pretty(range(agg.losses), n = nclass.FD(agg.losses), min.n = 1)
bwidth <- breaks[2]-breaks[1]
df <- data.frame(agg.losses)
ggplot(df,aes(agg.losses))+geom_histogram(binwidth=bwidth,fill="white",colour="deepskyblue4")
### Percentiles at which dist will be calculated
perc <- c(0.25, 0.5, 0.75, 0.8, 0.85, 0.9, 0.95, 0.98, 0.99, 0.995, 0.999)
### Show basic stats
gross <- c(quantile(agg.losses, probs = perc), mean(agg.losses), sd(agg.losses))
Percentile <- c(perc, "Mean", "SD")
results.table <- data.frame(Percentile, gross)
results.table
## Percentile gross
## 1 0.25 1875.651
## 2 0.5 4726.273
## 3 0.75 9782.279
## 4 0.8 11527.919
## 5 0.85 13894.378
## 6 0.9 17487.011
## 7 0.95 24421.672
## 8 0.98 35706.028
## 9 0.99 46308.625
## 10 0.995 58970.708
## 11 0.999 98352.501
## 12 Mean 7652.167
## 13 SD 10210.565