In sunhfsc/SC19052: Functions for collective risk model and risk measures

Overview

SC19052 is a simple R package developed to be used in collective risk model and risk measures.There are four functions in the R package,namely, ParaApp, Panjer.Poisson, SparseVec,and VaRplot,the three R functions can compute the cumulative loss for collective risk model,and a drawing function for risk measures can draw the curve of VaR and TVaR changing with quantile level.

Collective risk model calculation

In the collective risk model, it is assumed that the number of losses of an insurance policy during an insurance period is a random variable N, and the amount of the i-th loss is another random variable $X_{i}$, in which the loss amount of each time $X_{i}$ is independently and uniformly distributed. The cumulative loss of the collective risk model can be expressed as $$ S=X_{1}+X_{2}+\cdots+X_{N} $$ Then the exact distribution of cumulative loss is $$ \operatorname{Pr}(S=s)=\sum_{n=0}^{\infty} \operatorname{Pr}(S=s | N=n) \operatorname{Pr}(N=n) $$

The method of parametric approximation

In the collective risk model, the cumulative loss can be calculated by parametric approximation when N is enough large. The so-called parametric approximation is to approximate the cumulative loss distribution by a parameter distribution, mainly including the normal approximation, the translational gamma approximation, the normal power approximation and the Wilson-Hilfery approximation which as shown below,and $\mu$ , $\sigma$ and $\kappa$ represent the mean,standard deviation and skewness coefficients of the cumulative loss S, respectively. (1)Normal approximation $$ F_{s}(s)=\operatorname{Pr}(S \leqslant s)=\operatorname{Pr}\left(\frac{S-\mu}{\sigma} \leqslant \frac{s-\mu}{\sigma}\right) \approx \Phi\left(\frac{s-\mu}{\sigma}\right) $$

(2)Translational gamma approximation $$ \left{\begin{array}{l} {\alpha=\frac{4}{\kappa^{2}}} \ {\beta=\frac{2}{\kappa \sigma}} \ {x_{0}=\mu-\frac{2 \sigma}{\kappa}} \end{array}\right. $$ $$ \begin{aligned} F_{s}(s) &=\operatorname{Pr}(S \leqslant s) \approx \operatorname{Pr}\left(W+x_{0} \leqslant s\right) =\operatorname{Pr}\left(W \leqslant s-x_{0}\right) =F_{\mathrm{G}}\left(s-x_{0}\right) \end{aligned} $$

(3)Normal power approximation $$ F_{s}(s) \approx \Phi\left(-\frac{3}{\kappa}+\sqrt{\frac{9}{\kappa^{2}}+1+\frac{6}{\kappa} \frac{s-\mu}{\sigma}}\right) $$

(4)Wilson-Hilfery approximation $$ F_{s}(s) \approx \Phi\left[3\left(\frac{2}{\kappa}\right)^{2 / 3}\left(\frac{s-\mu}{\sigma}+\frac{2}{\kappa}\right)^{1 / 3}-\frac{6}{\kappa}+\frac{\kappa}{6}\right] $$

The function-"ParaApp" can approximate the cumulative loss distribution by choosing any of the four methods.
The source R code for ParaApp is as follows:

ParaApp <- function(s,mu,sigma,k,method){
  if(method=="normal"){
    Fs <- pnorm((s-mu)/sigma)
  }
  if(method=="Tgamma"){
    alpha <- 4/k^2
    beta <- 2/(k*sigma) 
    x0 <- mu-2*sigma/k
    Fs <- pgamma(s-x0,shape=alpha,rate=beta)
  }
  if(method=="normal power"){
    Fs <- pnorm(-3/k+sqrt(9/k^2+1+6/k*(s-mu)/sigma))
  }
  if(method=="Wilson Hilfery"){
    Fs <- pnorm(3*(2/k)^(2/3)*((s-mu)/sigma+2/k)^(1/3)-6/k+k/6)
  }
  return(Fs)
}

Example 1
Supposed that the number of losses N obeys the poisson distribution with a parameter of 3, and the amount of losses X obeys the exponential distribution with a ratio parameter of 0.01.
We can get the the mean, variance and skewness coefficients of the cumulative loss S: $$ \mu=300, \quad \sigma^{2}=60000, \quad \kappa=1.2247 $$

s <- c(0:1000)
mu <- 300
sigma <- sqrt(60000)
k <- 1.2247
Fs1 <- ParaApp(s,mu,sigma,k,method="normal")
Fs2 <- ParaApp(s,mu,sigma,k,method="Tgamma")
Fs3 <- ParaApp(s,mu,sigma,k,method="normal power")
Fs4 <- ParaApp(s,mu,sigma,k,method="Wilson Hilfery")
plot(s,Fs1,pch="",ylab="Fs")
lines(s,Fs1,col='red')
lines(s,Fs2,col='blue')
lines(s,Fs3,col='green')
lines(s,Fs4,col='yellow')

Panjer recursion

According to the theorem of Panjer recursion, when $N \sim Poisson(\lambda)$, Xi(i=1,2...N) has probability distribution function p(x),x=0,1,2...,then the probability distribution function of the cumulative loss S satisfies the following relation: $$ \begin{aligned} &f(0)=\mathrm{e}^{-\lambda(1-p(0))}\ &f(s)=\frac{1}{s} \sum_{h=1}^{s} \lambda h p(h) f(s-h) \end{aligned} $$

The function-"Panjer.Poisson" can be used to calculate the relation above.

Panjer.Poisson <- function(p,lambda){
  if(sum(p)>1||any(p<0)) stop("p parameter not a density")
  if(lambda*sum(p)>727) stop("Underflow")
  cumul <- f <- exp(-lambda*sum(p))
  r <- length(p)
  s <- 0
  repeat{
    s <- s+1
    m <- min(s,r)
    last <- lambda/s*sum(1:m*head(p,m)*rev(tail(f,m)))
    f <- c(f,last)
    cumul <- cumul+last
    if(cumul>0.99999999)
    break
  }
  return(list(f=f,s=s))
}

Example 2
Supposed that $N \sim Poisson(4)$,and $\operatorname{Pr}[X=1,2,3]=\frac{1}{4}, \frac{1}{2}, \frac{1}{4}$.

fs <- Panjer.Poisson(c(0.25,0.5,0.25),4)
s <- c(0:fs$s)
print(cbind(s,fs$f))

The method of sparse vector

The function-"SparseVec" also can compute the f(s)-the probability distribution function of the cumulative loss S.

SparseVec <- function(p,lambda){
  freq <- p*lambda
  if (any(freq<0)) stop("negative frequency")
  M <- length(freq)
  mu <- sum((1:M)*freq)
  sigma2 <- sum((1:M)^2*freq)
  MM <- ceiling(mu+10*sqrt(sigma2))+6
  fs <- dpois(0:(MM-1),freq[1])
  for(j in 2:M){
    MMM <- trunc((MM-1)/j)
    fj <- rep(0,MM)
    fj[(0:MMM)*j+1] <- dpois(0:MMM,freq[j])
    fs <- convolve(fs,rev(fj),type="o")
  }
  return(list(f=fs,s=length(fs)-1))
}

Example 3 Supposed that $N \sim Poisson(4)$,and $\operatorname{Pr}[X=1,2,3]=\frac{1}{4}, \frac{1}{2}, \frac{1}{4}$.

f <- SparseVec(c(0.25,0.5,0.25),4)
print(f)

Risk measures

In risk measurement, VaR(Value at Risk) is commonly used to describe the maximum possible loss at a certain confidence level. $$ \operatorname{VaR}{p}(X)=\inf {x \in R | F(x) \geqslant p}=\inf {x \in R | \bar{F}(x)<1-p} $$ TVaR is defined on the basis of VaR and refers to the expected value of the loss beyond VaR. The TVaR of the lost random variable X is defined as follows $$ \mathrm{TVaR}{p}(X)=E\left[X | X>\mathrm{VaR}{p}(X)\right] $$ The relation between them is $$ \operatorname{TVaR}{p}(X)=\operatorname{VaR}{p}(X)+\frac{P\left[X>\operatorname{VaR}{p}(X)\right]}{1-p} E\left[X-\operatorname{VaR}{p}(X) | X>\operatorname{VaR}{p}(X)\right] $$

The function-"VaRplot" can draw the curve of VaR and TVaR changing with quantile level.

VaRplot <- function(x){
p <- seq(0.01,0.99,0.001)
n <- length(p)

VaR=TVaR=NULL
for(j in 1:n){
  VaR[j] <- quantile(x,p[j])
  TVaR[j] <- VaR[j]+mean((x-VaR[j])*(x>VaR[j]))/(1-p[j])
}

plot(p,VaR,type="s",ylab="",lty=1)
lines(p,TVaR,type="s",lty=2)
legend("topleft",c("VaR","TVaR"),lty=c(1,2),bty="n")
}