The capital modeling approach often uses complex methods such as the Monte Carlo simulation to estimate with greater precision and confidence than a limited dataset might allow.

iihsData <- read.csv("C:/Users/joshu/Desktop/iihs.csv")

str(iihsData)
## 'data.frame':    110 obs. of  8 variables:
##  $ Vehicle         : Factor w/ 42 levels "Acura TSX","Buick Verano",..: 1 1 1 1 2 3 3 4 4 5 ...
##  $ YEAR            : int  2012 2011 2010 2009 2012 2012 2011 2010 2009 2010 ...
##  $ Collision.      : int  106 105 111 111 82 95 105 94 94 80 ...
##  $ Property.damage.: int  83 84 73 66 68 97 92 83 79 78 ...
##  $ Compre..hensive.: int  110 117 116 116 86 101 98 98 113 NA ...
##  $ Personal.injury.: int  94 93 89 88 114 122 115 111 103 NA ...
##  $ Medical.payment.: int  84 93 93 89 NA 136 119 112 104 NA ...
##  $ Bodily.injury.  : int  82 78 71 62 NA 123 114 105 93 NA ...
#Crashes Per Year
table(iihsData$YEAR)
## 
## 2009 2010 2011 2012 
##   28   29   27   26

Draw a Sample From a Poisson Distribution With a Mean Equal to the Average Number of Loss Events Per Year

table(iihsData$YEAR)*mean(iihsData$Collision)
## 
##     2009     2010     2011     2012 
## 2945.855 3051.064 2840.645 2735.436
yearMean <- mean(table(iihsData$YEAR))

#rpois is a Base R Function

rpois(10, yearMean) #10 Samples
##  [1] 33 27 30 19 28 31 26 29 35 16
rp1 <- rpois(1, yearMean) #1 Sample

rp1
## [1] 22

Determine the Average and Standard Deviation of the “Collision” Loss Amount

m <- mean(iihsData$Collision)

m
## [1] 105.2091
s <- sd(iihsData$Collision)

s
## [1] 14.83431
#Mean Log Normal Distribution
mu <- log(m^2/sqrt(s+m^2)) 

mu
## [1] 4.65528
sigma <- sqrt(log(1+s/m^2))

sigma
## [1] 0.03659612

Draw a Sample Value From a lognormal Distribution

rlnorm(10, mu, sigma)
##  [1] 107.96242 108.61045  99.72251 110.48398 110.20688 107.92588 108.04458
##  [8] 109.26818 101.26600 104.54957
n1 <- rlnorm(1, mu, sigma)

n1
## [1] 108.2607

Multiply the random value drawn from the poisson distribution by the random value drawn from the lognormal distribution.

f1 <- rp1 * n1

Repeat the process four more times.

rp2 <- rpois(1, yearMean) 
rp3 <- rpois(1, yearMean) 
rp4 <- rpois(1, yearMean) 
rp5 <- rpois(1, yearMean) 

n2 <- rlnorm(1, mu, sigma)
n3 <- rlnorm(1, mu, sigma)
n4 <- rlnorm(1, mu, sigma)
n5 <- rlnorm(1, mu, sigma)

f2 <- rp2 * n2
f3 <- rp3 * n3
f4 <- rp4 * n4
f5 <- rp5 * n5

samples <- c(f1,f2,f3,f4,f5)

plot(samples)

Let’s try a loop.

#create a vector to store rp objects
rpObjects = c(rp1, rp2, rp3, rp4, rp5)

#loop to add 500 more members to rpObjects
for(i in 6:505)
  rpObjects[i] = rpois(1, yearMean)

str(rpObjects)
##  int [1:505] 22 25 29 23 28 30 31 27 30 42 ...
#create a vector to store n objects
nObjects = c(n1, n2, n3, n4, n5)

for(i in 6:505)
  nObjects[i] = rlnorm(1, mu, sigma)

str(nObjects)
##  num [1:505] 108 113 102 113 107 ...