Chapter 5 Simulation

This file contains illustrative R code for computing important count distributions. When reviewing this code, you should open an R session, copy-and-paste the code, and see it perform. Then, you will be able to change parameters, look up commands, and so forth, as you go.

5.1 Simulation - Inversion method

This section shows how to use the inversion method to simulate claims from a gamma distribution. The results below are summary statistics from the simulated data.

##Simulation-gamma
library(moments)
set.seed(2) #set seed to reproduce work 
nTot<-20000  #number of simulations
alpha<-2
theta<-100
         
Losses<-rgamma(nTot,alpha,scale = theta)  
summary(Losses)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##    0.0921   96.3300  167.8000  200.2000  268.2000 1110.0000

k<-0.95
PercentileLoss<-quantile(Losses,k)     # Kth percentile of Losses 
PercentileLoss

##      95% 
## 473.8218

#######################################
#OR you can use this method to simulate losses 
#Fx<-runif(nTot)
#Losses<-qgamma(Fx,alpha,scale=theta)
#######################################
# For the Pareto Distribution, use
#library(VGAM)
# nTot<-10000  #number of simulations
# alpha<-3
# theta<-25000
# Losses<-rparetoII(nTot,scale=theta,shape=alpha)
# rparetoII(nTot,scale=theta,shape=alpha)

A few quick notes on these commands:

The rgamma() function randomly generates data from the Gamma distribution. In this illustration the data was generated from a gamma distribution with parameters shape=alpha=2 and scale=theta=100.
The quantile() function provides sample quantiles corresponding to the given probabilities. Here we wanted the simulated loss data corresponding to the 95th percentile.

5.2 Comparing moments from the simulated data to theoretical moments

library(pander)
###Raw moments for k=0.5
# Theoretical 
k<-0.5
T0.5<-round(((theta^k)*gamma(alpha+k))/gamma(alpha),2)

#Simulated data raw moments
S0.5<-round(moment(Losses, order = k, central = FALSE),2)


###Raw moments for k=1
# Theoretical 
k<-1
T1<-((theta^k)*gamma(alpha+k))/gamma(alpha)

#Simulated data raw moments
S1<-round(moment(Losses, order = k, central = FALSE),2)

###Raw moments for k=2
# Theoretical 
k<-2
T2<-((theta^k)*gamma(alpha+k))/gamma(alpha)

#Simulated data raw moments
S2<-round(moment(Losses, order = k, central = FALSE),2)

###Raw moments for k=3
# Theoretical 
k<-3
T3<-((theta^k)*gamma(alpha+k))/gamma(alpha)

#Simulated data raw moments
S3<-round(moment(Losses, order = k, central = FALSE),2)

###Raw moments for k=3
# Theoretical 
k<-4
T4<-((theta^k)*gamma(alpha+k))/gamma(alpha)

#Simulated data raw moments
S4<-round(moment(Losses, order = k, central = FALSE),2)

pander(rbind(c("k",0.5,1,2,3,4),c("Theoritical",T0.5,T1,T2,T3,T4),c("Simulated",S0.5,S1,S2,S3,S4)))

k	0.5	1	2	3	4
Theoritical	13.29	200	60000	2.4e+07	1.2e+10
Simulated	13.3	200.17	60069.73	23993892.53	11897735665.89

A few quick notes on these commands:

The momemt() function computes all the sample moments of the chosen type up to a given order. In this illustration, we wanted the raw sample moments of the kth order
The round() function was used to round values.