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] 26 26 27 22 25 30 36 24 23 25
rp1 <- rpois(1, yearMean) #1 Sample

rp1
## [1] 25

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] 104.8696 102.0647 102.2849 106.7471 106.2142 111.9814 108.2703
##  [8] 100.9169 105.3986 107.2610
n1 <- rlnorm(1, mu, sigma)

n1
## [1] 107.8723

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 95 more members to rpObjects
for(i in 6:100)
  rpObjects[i] = rpois(1, yearMean)

str(rpObjects)
##  int [1:100] 25 25 35 26 22 26 32 28 35 25 ...
rpMean <- mean(rpObjects)

rpMean
## [1] 28.3
rpSd <- sd(rpObjects)

rpSd
## [1] 5.40015
#create a vector to store n objects
nObjects = c(n1, n2, n3, n4, n5)

#similar loop
for(i in 6:100)
  nObjects[i] = rlnorm(1, mu, sigma)

str(nObjects)
##  num [1:100] 107.9 103.9 107 108.2 99.4 ...
nMean <- mean(nObjects)

nMean
## [1] 105.8502
nSd <- sd(nObjects)

nSd
## [1] 3.867755
#create a vector to hold f objects
fObjects = c(f1, f2, f3, f4, f5)

#similar loop
for(i in 6:100)
  fObjects[i] = rpObjects[i] * nObjects[i]
  
str(fObjects)
##  num [1:100] 2697 2598 3747 2813 2188 ...
hist(fObjects)

fMean <- mean(fObjects)

fMean
## [1] 2997.127
fSd <- sd(fObjects)

fSd
## [1] 589.3491

Simulations for # tries for normal distribution and CI with these functions

rpObjects2 = c(rp1, rp2, rp3, rp4, rp5)
nObjects2 = c(n1, n2, n3, n4, n5)
fObjects2 = c(f1, f2, f3, f4, f5)

#try by 100 million, looks pretty normal

for(i in 6:100000)
  rpObjects2[i] = rpois(1, yearMean)

for(i in 6:100000)
  nObjects2[i] = rlnorm(1, mu, sigma)

for(i in 6:100000)
  fObjects2[i] = rpObjects2[i]*nObjects2[i]

str(fObjects2)
##  num [1:100000] 2697 2598 3747 2813 2188 ...
plot(fObjects2)

hist(fObjects2)

mm <- mean(fObjects2)

mm
## [1] 2890.904
ss <- sd(fObjects2)

ss
## [1] 562.2601