# Pricing Zero-Coupon Catastrophe Bonds Using EVT with Doubly Stochastic Poisson Arrivals

Authors: Zonggang Ma; Chaoqun Ma; Shisong Xiao
Journal: Discrete Dynamics in Nature and Society
 (2017)

---

## Abstract

The frequency and severity of climate abnormal change displays an irregular upward cycle as global warming intensifies. Therefore, this paper employs a doubly stochastic Poisson process with Black Derman Toy (BDT) intensity to describe the catastrophic characteristics. By using the Property Claim Services (PCS) loss index data from 2001 to 2010 provided by the US Insurance Services Office (ISO), the empirical result reveals that the BDT arrival rate process is superior to the nonhomogeneous Poisson and lognormal intensity process due to its smaller RMSE, MAE, MRPE, and U and larger E and d. Secondly, to depict extreme features of catastrophic risks, this paper adopts the Peak Over Threshold (POT) in extreme value theory (EVT) to characterize the tail characteristics of catastrophic loss distribution. And then the loss distribution is analyzed and assessed using a quantile-quantile (QQ) plot to visually check whether the PCS index observations meet the generalized Pareto distribution (GPD) assumption. Furthermore, this paper derives a pricing formula for zero-coupon catastrophe bonds with a stochastic interest rate environment and aggregate losses generated by a compound doubly stochastic Poisson process under the forward measure. Finally, simulation results verify pricing model predictions and show how catastrophic risks and interest rate risk affect the prices of zero-coupon catastrophe bonds.


---

## Body

## 1. Introduction

The Swiss Re Sigma world insurance database shows that the significant increases of “what" over the past several years are not due to an upward trend in the frequency and severity of natural catastrophic risk events but rather to increasing concentrations of population and property values in catastrophe prone areas since the late 1980s, as revealed in Figures1, 2, and 3. The Property Claim Services (PCS) data present an irregular upward cycle in Figures 4 and 5. The Intergovernmental Panel on Climate Change (IPCC) Fifth Assessment Report (2013) also forecasts that the frequency of global extreme disasters will continue to increase in the 21st century, and this seriously influences the sustainable development of the economy and the stability of society [1]. And faced with severe global climate changes, insurers might not have sufficient reinsurance purchases, or reinsurance providers might not have made sufficient capital to satisfy their existing obligations. In any catastrophic risk event (such as Hurricane Andrew in 1992 and Hurricanes Sandy in 2012), after a catastrophic risk loss, the capital of reinsurer would be seriously weakened; that is, new company formation through initial public offerings was insufficient to enable the insurance market to recover to previous levels of reinsurance capacity [2, 3].Figure 1
Number of global natural catastrophes from 1970 to 2016.Figure 2
Adjusted total losses from natural catastrophes worldwide during 1970–2016 (USD billion, 2016 prices).Figure 3
Adjusted global insured claims from natural catastrophes during 1970–2016 (USD billion, 2016 prices).Figure 4
Number of natural CATs in the United States from 2001 to 2010.Figure 5
The PCS catastrophe loss data in the United States from 2001 to 2010 (USD billion, 2010 prices).Alternative risk transfer (also called ART) has become increasingly more important over the past two decades, especially due to an increasing risk of extreme losses caused by irregular climate change and ascribed to a limited capacity of traditional reinsurance markets [2]. In response to a capacity shortage in reinsurance markets, ART intends to provide additional (re)insurance coverage by transferring (re)insurance risks to the powerful capital market, which offers considerably higher capacities and can help meet the demand. The most prominent type of ART is catastrophe risk (CAT) bonds, which are fully collateralized financial instrument that pays off on the occurrence of a preassigned catastrophic risk event [4]. The first ever catastrophe bond was a USD 85 million issue by Hanover Re in 1994 [5]. The market for CAT bonds has developed rapidly since their introduction about 20 years ago. Heading into 2016, the catastrophe bonds market continued to provide an increasingly attractive and valuable supplement to sponsor risk transfer programs; the Artemis Deal Directory (http://www.artemis.bm/deal_directory/) shows that USD 6.293 billion in risk capital was issued through 37 transactions in 2016, and the volume of CAT bonds principal outstanding rose to USD 25.752 billion due on March 2017.Accurate pricing is the key point in the CAT bonds issuing and trading. Recently, some studies have mainly focused on catastrophe bonds pricing models involved in a compound doubly stochastic Poisson loss process to depict the dynamic aggregate losses accurately [6, 7]. Nevertheless, the assumption that resulting losses occur in terms of the Poisson loss process has a constant arrival rate. As we all know, for unanticipated catastrophic risk events, the frequency and severity of climate abnormal change follow an irregular upward cycle as the global warms. Therefore, it seems insufficient to assume that catastrophic risk arrival process follows a Poisson process with a deterministic intensity function. This paper utilizes a doubly stochastic Poisson process with stochastic intensity to model the catastrophic risk arrival process. Furthermore, this paper employs the catastrophe loss index provided by Property Claim Services (PCS) from 2001 to 2010 to test the quality of the fitting between the stochastic intensity and the deterministic intensity. Empirical result shows that the doubly stochastic Poisson process is appropriate to model the arrival process for catastrophic risk events. This study’s primary contribution is to provide a semianalytical solution for evaluating zero-coupon CAT bonds with a doubly stochastic compound Poisson process pricing framework under the forward measure.The paper is organized as follows. Section2 describes the model assumptions and illustrates CAT claim models. In Section 3 we show the forward-neutral measure and derive the pricing formula for zero-coupon CAT bonds. Section 4 is devoted to empirical experiments. Section 5 presents a numerical analysis. Finally, conclusions are presented in Section 6.


## 2. Valuation Methodology and the Modeling Assumptions

### 2.1. Setup and Notation

For ease of exposition, we first introduce the basic setup, which will be generalized in the following sections. Let{VCAT(t):t∈[0,T]} denote the CAT bond price process, which is modeled by many factors: type of region, kind of loss event, sort of insured property, interest rate uncertainty, and so on. Let {Lt:t∈[0,T]} denote the aggregate loss process of the insured up to time t. {Nt:t∈[0,T]} is a doubly stochastic Poisson point process with an intensity parameter λt. {Xj:j⩾1} is a sequence of identically and independently distributed (i.i.d.) random variables and accounts for the sizes of CAT loss.In addition, let{M(mt):t∈[0,T]} stand for a doubly stochastic Poisson process with stochastic intensity mt and represent the CAT number, and let {r(t):t∈[0,T]} denote the risk-free short rate process (or the force of interest). {Wm(t):t∈[0,T]} and {Wr(t):t∈[0,T]} denote the Wiener process and account for two sources of randomness of the intensity rate of CAT occurrence and the uncertainty of interest rates, respectively. The interest rate process and the loss process are assumed to be stochastically independent. Let (Ω,F,p) define a probability space, where Ω is the set of states of the world, the filtration F is σ-algebra of subsets of Ω and is defined by {Ft:t∈[0,T]}, where Ft is σ(Wr(s),Wm(s),Xj:0≤s≤t,j≥1) and P is a probability measure on Ft.LetD(t,T) denote the price at time t of a (default-free) discount bond that pays USD 1 at maturity date T, which is given by(1)Dt,T=EQe-∫tTrsds∣Ft.Then, the spot forward rate at time t for future date T is defined by(2)ft,T=-∂log⁡Dt,T∂T.


### 2.2. Valuation Theory

In a risk-free arbitrage financial market, the priceV(t,Xt) of the contingent claim ψ(XT) at time t based on an underlying asset X at maturity T is given by(3)Vt,Xt=EQe-∫tTrsdsψXT∣Ft,where EQ denotes the conditional expectation and ψ(XT) is the payoff of the contingent claim. For more information about the fundamental theorem of asset pricing, please refer to Harrison and Kreps [8] and Harrison and Pliska [9].An indispensable point is the choice of the numeraire when determining the forward-neutral measureQT. That is, the common unit on the basis of which asset prices is expressed. Any asset price can be selected as a numeraire, as long as it has a strictly positive value in any state of the world. A very popular choice is to adopt the T-maturity bond D(t,T) as a numeraire for contingent claims that have a payoff ψ(XT) at time T. By applying the Change of Numeraire Theorem, the prices V(t,Xt)/D(t,T) are martingales for t<T under the measure QT. Hence, by applying the definition of a martingale, we obtain(4)EQTVT,XTDT,T∣Ft=Vt,XtDt,T.However, at maturity T, the price of the discount bond D(T,T)≡1 and the price of the contingent claim are given by the payoff V(T,XT)=ψ(XT). So, (4) reduces to(5)Vt,Xt=Dt,TEQTVT,XT∣Ft=Dt,TEQTψXT∣Ft.Taking derivatives of (1) with respect to T yields(6)-∂Dt,T∂T=EQe-∫tTrsdsrT∣Ft=Dt,TEQe-∫tTrsdsDt,TrT∣Ft=Dt,TEQTrT∣Ft,where we have used the Change of Numeraire Theorem in the last step with Randon-Nikodym derivative(7)dQTdQ=e-∫tTrsdsDt,T.Then we can simplify (2) to(8)ft,T=EQTrT∣Ft.Therefore, the spot forward ratef is a martingale under the measure QT.


### 2.3. Interest Rate Dynamics

There is no doubt that the interest rate risk is one of the important risk sources for both the insurer and reinsurer. To assess the effect of interest rate risk on the catastrophe bond prices, most previous studies use a constant interest rate or single-factor equilibrium model, such as Burnecki and Kukla (2003) [6] and Härdle and Cabrera (2010) [7]. However, a disadvantage of such models is that they give a set of theoretical prices for bonds that will not normally precisely match the actual prices that we observe in the markets. Therefore, this paper assumes that the (re)insurer operates in an environment where interest rates are stochastic and follow the Hull-White model [10]. The spot interest rate process is denoted by the following stochastic differential equation (SDE): (9)drt=θt-κrtdt+σrdWrt,where κ and σr are positive constants and θ(t) is a deterministic function of time. The parameter θ(t) here can be analytically computed as (10)θt=∂f0,t∂t+κf0,t+σr22κ1-e-2κt.Letr(t)=y(t)+β(t) and β(t)=e-κt(r(0)+∫0teκuθ∗(u)du). Then, the mean and variance of y(T) are given by (11)meanyT=e-κT-tyt,VaryT=σr22κ1-e-2κT-t,for 0≤t≤T.


### 2.4. Aggregate Claim Dynamics

The aggregate claim has a compound distribution with two main components: the frequency and severity of catastrophic risk events [11].


#### 2.4.1. The Claim Arrival Process

As we all know, the pure Poisson process has been frequently used as a claim arrival process, which generally assumes that the number of claims observed until timet is a (non)homogeneous Poisson process, denoted by Mt, with the intensity mt of the counting process varying with time t (e.g., see [12, 13]). However, for unanticipated catastrophic events, the assumption that resulting claims occur in terms of a pure Poisson process is insufficient because it has deterministic intensity [14]. Research has shown that the doubly stochastic Poisson process is more appropriately used as a claim arrival process as it can explain the assumption that catastrophic risk events occur periodically [15].Let{Mt} be a point process adapted to Ft, and let mt be a nonnegative process. Suppose that mt is Ft- measurable, t≥0, and that (12)∫0tmsds<∞almost surely (noexplosions).If for all 0≤t1≤t2 and u∈R(13)EeiuMt2-Mt1∣Ft1=exp⁡eiu-1∫t1t2msds,then Mt is called Ft-doubly stochastic Poisson process with intensity mt.Equation (13) yields, for all t1≤s≤t2 and k≥0, (14)PMt2-Mt1=k∣ms;t1≤s≤t2=e-∫t1t2msds∫t1t2msdskk!,where k! denotes a factorial.Thus, the probability of no occurrences within the interval[t1,t2] for the process Mt with intensity mt>0 is given by (15)PMt2-Mt1=0=Eexp⁡-∫t1t2msds.Accordingly, to capture more features of the catastrophic process, we employ a doubly stochastic Poisson process to describe the catastrophic events arrival process, whose intensity follows the Black, Derman, and Toy (BDT) model [16], which is given by the following SDE: (16)dln⁡mt=αtdt+νdWmt,where α(t) is a deterministic function of t and ν is a positive constant.In addition, to describe a complex periodic characteristic of catastrophic risk events, this paper takes hereα(t)=a+bsin⁡(π(t+c)), where a≥0,b and c are constants.By integrating both sides of (16), this paper obtains the exact solution of the instantaneous intensity: (17)ln⁡mt=ln⁡m0+∫0tαsds+∫0tνdWms=ln⁡m0+∫0tαsds+νWms.Then, (18)mt=m0exp⁡∫0tαsds+νWmt,or for all 0≤u≤t(19)mt=muexp⁡∫utαsds+νWmt-Wmu,where αt=a+bsin⁡(π(t+c)),a≥0, andb and c are constants. Figure 6 presents an example of a trajectory followed by the BDT intensity process. It can be observed that these curves exhibit certain periodic and upward trends.Figure 6
Example of path for the stochastic intensitymt. The graph corresponds to the parameter values m0=20, a=0.05, b=-0.3, and c=-0.5.


#### 2.4.2. Modeling Loss Severities

Extreme value theory (EVT) provides a simple technique for evaluating probabilities of future extreme levels of a process based on historical observations. Generally, being a large claim means exceeding some prespecified threshold value. The generalized Pareto distribution (GPD) can be used to set the value of this threshold (see for a detailed account in the book by Embrechts et al. [17]). Suppose insured claims are denoted by random variables X over a high threshold u with distribution function F. Then the distribution is defined by(20)Fuy=pX≤u+y∣X>u=Fy+u-Fu1-Fu=Fx-Fu1-Fu,which denotes the probability that the value of x exceeds the threshold u by at most an amount y given that x exceeds the threshold u. Balkema and de Haan [18] and Pickands III [19] showed that, for sufficiently large threshold u, the distribution function of the excess can be approximated by the GPD such that, as the threshold u gets high, the excess distribution Fu(y) converges to the GPD, which is(21)Gξ,σXx=1-1+ξxσX-1/ξifξ≠0,1-exp⁡-xσXifξ=0,where σX>0, and the range for x is x≥0, when ξ≥0, and x∈[0,-σX/ξ] when ξ<0. The parameters σ and ξ are expressed as the scale and shape parameters, respectively.Suppose thatX={X1,X2,…,Xnu} is a random sampling of the GPD, where nu denote the sample exceeding a threshold level u. For each of these exceedances we calculate the amount yi=Xi-u(i=1,2,…,nu) of the excess loss. These generalized Pareto distribution parameters may be evaluated by maximum likelihood estimation (MLE) using the following stages:(1)
Make log-likelihood function forξ≠0(22)lθ∣y=ln⁡Lθ∣y=-nuln⁡σX-∑i=1nu1+1ξln⁡1+ξyiσX,whereθ=(ξ,σX), σX>0, 1+ξyi/σX>0.(2)
Solve the maximization of (22) by using numerical optimization algorithms. Through differentiating (22) with respect to ξ and σX and then equating to zero, this paper gets the following equations: (23)∂lθ∣y∂ξ=1ξ2∑i=1nuln⁡1+ξyiσX+ξyi-∑i=1nuyiσX+ξyi=0,∂lθ∣y∂σX=-nuσX+1ξ+1∑i=1nuξyiσX2+σXξyi=0.


#### 2.4.3. Doubly Stochastic Compound Poisson Process

The risk faced by the (re)insurance industry is intrinsic in their exposure to aggregate claims. This paper assumes that the process of aggregate claims is a compound doubly stochastic Poisson process. Then we need to make some assumptions about the aggregate claims process.(i)
There exists a doubly stochastic Poisson processMs(s∈[0,T]) representing potentially catastrophic events. The Poisson point arrival process is assumed to be a predictable bounded process ms. This paper denotes the moments of these potentially catastrophic risk as 0≤t1≤⋯≤⋯≤tj≤⋯≤T.(ii)
The insured claims incurred by each catastrophic risk in the flow{tj}j=1,… are assumed to be independently and identically distributed random variables {Xj}j=1,… with a distribution function F(x)=P{Xj<x} and a density function f(x).Therefore, a left-continuous and predictable aggregate claims process can be defined as(24)Lt=∑tj<tXj=∑j=1MtXj,where Mt and Xj are assumed to be independent.In addition, there is also another assumption:(iii)
The threshold level is the timeτ when the aggregate claims Lt exceeds the threshold value D; then, τ=inf{t:Lt≥D}. Now defining a new process Nt=1{Lt≥D}, Baryshnikov et al. [20] showed that it is a doubly stochastic Poisson process with stochastic intensity: (25)λs=ms1-FD-Ls1Ls<D.However, in that paper, the proof ofλt is rather loose. So, this paper will give a detailed proof of the result in Appendix A.


## 2.1. Setup and Notation

For ease of exposition, we first introduce the basic setup, which will be generalized in the following sections. Let{VCAT(t):t∈[0,T]} denote the CAT bond price process, which is modeled by many factors: type of region, kind of loss event, sort of insured property, interest rate uncertainty, and so on. Let {Lt:t∈[0,T]} denote the aggregate loss process of the insured up to time t. {Nt:t∈[0,T]} is a doubly stochastic Poisson point process with an intensity parameter λt. {Xj:j⩾1} is a sequence of identically and independently distributed (i.i.d.) random variables and accounts for the sizes of CAT loss.In addition, let{M(mt):t∈[0,T]} stand for a doubly stochastic Poisson process with stochastic intensity mt and represent the CAT number, and let {r(t):t∈[0,T]} denote the risk-free short rate process (or the force of interest). {Wm(t):t∈[0,T]} and {Wr(t):t∈[0,T]} denote the Wiener process and account for two sources of randomness of the intensity rate of CAT occurrence and the uncertainty of interest rates, respectively. The interest rate process and the loss process are assumed to be stochastically independent. Let (Ω,F,p) define a probability space, where Ω is the set of states of the world, the filtration F is σ-algebra of subsets of Ω and is defined by {Ft:t∈[0,T]}, where Ft is σ(Wr(s),Wm(s),Xj:0≤s≤t,j≥1) and P is a probability measure on Ft.LetD(t,T) denote the price at time t of a (default-free) discount bond that pays USD 1 at maturity date T, which is given by(1)Dt,T=EQe-∫tTrsds∣Ft.Then, the spot forward rate at time t for future date T is defined by(2)ft,T=-∂log⁡Dt,T∂T.


## 2.2. Valuation Theory

In a risk-free arbitrage financial market, the priceV(t,Xt) of the contingent claim ψ(XT) at time t based on an underlying asset X at maturity T is given by(3)Vt,Xt=EQe-∫tTrsdsψXT∣Ft,where EQ denotes the conditional expectation and ψ(XT) is the payoff of the contingent claim. For more information about the fundamental theorem of asset pricing, please refer to Harrison and Kreps [8] and Harrison and Pliska [9].An indispensable point is the choice of the numeraire when determining the forward-neutral measureQT. That is, the common unit on the basis of which asset prices is expressed. Any asset price can be selected as a numeraire, as long as it has a strictly positive value in any state of the world. A very popular choice is to adopt the T-maturity bond D(t,T) as a numeraire for contingent claims that have a payoff ψ(XT) at time T. By applying the Change of Numeraire Theorem, the prices V(t,Xt)/D(t,T) are martingales for t<T under the measure QT. Hence, by applying the definition of a martingale, we obtain(4)EQTVT,XTDT,T∣Ft=Vt,XtDt,T.However, at maturity T, the price of the discount bond D(T,T)≡1 and the price of the contingent claim are given by the payoff V(T,XT)=ψ(XT). So, (4) reduces to(5)Vt,Xt=Dt,TEQTVT,XT∣Ft=Dt,TEQTψXT∣Ft.Taking derivatives of (1) with respect to T yields(6)-∂Dt,T∂T=EQe-∫tTrsdsrT∣Ft=Dt,TEQe-∫tTrsdsDt,TrT∣Ft=Dt,TEQTrT∣Ft,where we have used the Change of Numeraire Theorem in the last step with Randon-Nikodym derivative(7)dQTdQ=e-∫tTrsdsDt,T.Then we can simplify (2) to(8)ft,T=EQTrT∣Ft.Therefore, the spot forward ratef is a martingale under the measure QT.


## 2.3. Interest Rate Dynamics

There is no doubt that the interest rate risk is one of the important risk sources for both the insurer and reinsurer. To assess the effect of interest rate risk on the catastrophe bond prices, most previous studies use a constant interest rate or single-factor equilibrium model, such as Burnecki and Kukla (2003) [6] and Härdle and Cabrera (2010) [7]. However, a disadvantage of such models is that they give a set of theoretical prices for bonds that will not normally precisely match the actual prices that we observe in the markets. Therefore, this paper assumes that the (re)insurer operates in an environment where interest rates are stochastic and follow the Hull-White model [10]. The spot interest rate process is denoted by the following stochastic differential equation (SDE): (9)drt=θt-κrtdt+σrdWrt,where κ and σr are positive constants and θ(t) is a deterministic function of time. The parameter θ(t) here can be analytically computed as (10)θt=∂f0,t∂t+κf0,t+σr22κ1-e-2κt.Letr(t)=y(t)+β(t) and β(t)=e-κt(r(0)+∫0teκuθ∗(u)du). Then, the mean and variance of y(T) are given by (11)meanyT=e-κT-tyt,VaryT=σr22κ1-e-2κT-t,for 0≤t≤T.


## 2.4. Aggregate Claim Dynamics

The aggregate claim has a compound distribution with two main components: the frequency and severity of catastrophic risk events [11].


### 2.4.1. The Claim Arrival Process

As we all know, the pure Poisson process has been frequently used as a claim arrival process, which generally assumes that the number of claims observed until timet is a (non)homogeneous Poisson process, denoted by Mt, with the intensity mt of the counting process varying with time t (e.g., see [12, 13]). However, for unanticipated catastrophic events, the assumption that resulting claims occur in terms of a pure Poisson process is insufficient because it has deterministic intensity [14]. Research has shown that the doubly stochastic Poisson process is more appropriately used as a claim arrival process as it can explain the assumption that catastrophic risk events occur periodically [15].Let{Mt} be a point process adapted to Ft, and let mt be a nonnegative process. Suppose that mt is Ft- measurable, t≥0, and that (12)∫0tmsds<∞almost surely (noexplosions).If for all 0≤t1≤t2 and u∈R(13)EeiuMt2-Mt1∣Ft1=exp⁡eiu-1∫t1t2msds,then Mt is called Ft-doubly stochastic Poisson process with intensity mt.Equation (13) yields, for all t1≤s≤t2 and k≥0, (14)PMt2-Mt1=k∣ms;t1≤s≤t2=e-∫t1t2msds∫t1t2msdskk!,where k! denotes a factorial.Thus, the probability of no occurrences within the interval[t1,t2] for the process Mt with intensity mt>0 is given by (15)PMt2-Mt1=0=Eexp⁡-∫t1t2msds.Accordingly, to capture more features of the catastrophic process, we employ a doubly stochastic Poisson process to describe the catastrophic events arrival process, whose intensity follows the Black, Derman, and Toy (BDT) model [16], which is given by the following SDE: (16)dln⁡mt=αtdt+νdWmt,where α(t) is a deterministic function of t and ν is a positive constant.In addition, to describe a complex periodic characteristic of catastrophic risk events, this paper takes hereα(t)=a+bsin⁡(π(t+c)), where a≥0,b and c are constants.By integrating both sides of (16), this paper obtains the exact solution of the instantaneous intensity: (17)ln⁡mt=ln⁡m0+∫0tαsds+∫0tνdWms=ln⁡m0+∫0tαsds+νWms.Then, (18)mt=m0exp⁡∫0tαsds+νWmt,or for all 0≤u≤t(19)mt=muexp⁡∫utαsds+νWmt-Wmu,where αt=a+bsin⁡(π(t+c)),a≥0, andb and c are constants. Figure 6 presents an example of a trajectory followed by the BDT intensity process. It can be observed that these curves exhibit certain periodic and upward trends.Figure 6
Example of path for the stochastic intensitymt. The graph corresponds to the parameter values m0=20, a=0.05, b=-0.3, and c=-0.5.


### 2.4.2. Modeling Loss Severities

Extreme value theory (EVT) provides a simple technique for evaluating probabilities of future extreme levels of a process based on historical observations. Generally, being a large claim means exceeding some prespecified threshold value. The generalized Pareto distribution (GPD) can be used to set the value of this threshold (see for a detailed account in the book by Embrechts et al. [17]). Suppose insured claims are denoted by random variables X over a high threshold u with distribution function F. Then the distribution is defined by(20)Fuy=pX≤u+y∣X>u=Fy+u-Fu1-Fu=Fx-Fu1-Fu,which denotes the probability that the value of x exceeds the threshold u by at most an amount y given that x exceeds the threshold u. Balkema and de Haan [18] and Pickands III [19] showed that, for sufficiently large threshold u, the distribution function of the excess can be approximated by the GPD such that, as the threshold u gets high, the excess distribution Fu(y) converges to the GPD, which is(21)Gξ,σXx=1-1+ξxσX-1/ξifξ≠0,1-exp⁡-xσXifξ=0,where σX>0, and the range for x is x≥0, when ξ≥0, and x∈[0,-σX/ξ] when ξ<0. The parameters σ and ξ are expressed as the scale and shape parameters, respectively.Suppose thatX={X1,X2,…,Xnu} is a random sampling of the GPD, where nu denote the sample exceeding a threshold level u. For each of these exceedances we calculate the amount yi=Xi-u(i=1,2,…,nu) of the excess loss. These generalized Pareto distribution parameters may be evaluated by maximum likelihood estimation (MLE) using the following stages:(1)
Make log-likelihood function forξ≠0(22)lθ∣y=ln⁡Lθ∣y=-nuln⁡σX-∑i=1nu1+1ξln⁡1+ξyiσX,whereθ=(ξ,σX), σX>0, 1+ξyi/σX>0.(2)
Solve the maximization of (22) by using numerical optimization algorithms. Through differentiating (22) with respect to ξ and σX and then equating to zero, this paper gets the following equations: (23)∂lθ∣y∂ξ=1ξ2∑i=1nuln⁡1+ξyiσX+ξyi-∑i=1nuyiσX+ξyi=0,∂lθ∣y∂σX=-nuσX+1ξ+1∑i=1nuξyiσX2+σXξyi=0.


### 2.4.3. Doubly Stochastic Compound Poisson Process

The risk faced by the (re)insurance industry is intrinsic in their exposure to aggregate claims. This paper assumes that the process of aggregate claims is a compound doubly stochastic Poisson process. Then we need to make some assumptions about the aggregate claims process.(i)
There exists a doubly stochastic Poisson processMs(s∈[0,T]) representing potentially catastrophic events. The Poisson point arrival process is assumed to be a predictable bounded process ms. This paper denotes the moments of these potentially catastrophic risk as 0≤t1≤⋯≤⋯≤tj≤⋯≤T.(ii)
The insured claims incurred by each catastrophic risk in the flow{tj}j=1,… are assumed to be independently and identically distributed random variables {Xj}j=1,… with a distribution function F(x)=P{Xj<x} and a density function f(x).Therefore, a left-continuous and predictable aggregate claims process can be defined as(24)Lt=∑tj<tXj=∑j=1MtXj,where Mt and Xj are assumed to be independent.In addition, there is also another assumption:(iii)
The threshold level is the timeτ when the aggregate claims Lt exceeds the threshold value D; then, τ=inf{t:Lt≥D}. Now defining a new process Nt=1{Lt≥D}, Baryshnikov et al. [20] showed that it is a doubly stochastic Poisson process with stochastic intensity: (25)λs=ms1-FD-Ls1Ls<D.However, in that paper, the proof ofλt is rather loose. So, this paper will give a detailed proof of the result in Appendix A.


## 2.4.1. The Claim Arrival Process

As we all know, the pure Poisson process has been frequently used as a claim arrival process, which generally assumes that the number of claims observed until timet is a (non)homogeneous Poisson process, denoted by Mt, with the intensity mt of the counting process varying with time t (e.g., see [12, 13]). However, for unanticipated catastrophic events, the assumption that resulting claims occur in terms of a pure Poisson process is insufficient because it has deterministic intensity [14]. Research has shown that the doubly stochastic Poisson process is more appropriately used as a claim arrival process as it can explain the assumption that catastrophic risk events occur periodically [15].Let{Mt} be a point process adapted to Ft, and let mt be a nonnegative process. Suppose that mt is Ft- measurable, t≥0, and that (12)∫0tmsds<∞almost surely (noexplosions).If for all 0≤t1≤t2 and u∈R(13)EeiuMt2-Mt1∣Ft1=exp⁡eiu-1∫t1t2msds,then Mt is called Ft-doubly stochastic Poisson process with intensity mt.Equation (13) yields, for all t1≤s≤t2 and k≥0, (14)PMt2-Mt1=k∣ms;t1≤s≤t2=e-∫t1t2msds∫t1t2msdskk!,where k! denotes a factorial.Thus, the probability of no occurrences within the interval[t1,t2] for the process Mt with intensity mt>0 is given by (15)PMt2-Mt1=0=Eexp⁡-∫t1t2msds.Accordingly, to capture more features of the catastrophic process, we employ a doubly stochastic Poisson process to describe the catastrophic events arrival process, whose intensity follows the Black, Derman, and Toy (BDT) model [16], which is given by the following SDE: (16)dln⁡mt=αtdt+νdWmt,where α(t) is a deterministic function of t and ν is a positive constant.In addition, to describe a complex periodic characteristic of catastrophic risk events, this paper takes hereα(t)=a+bsin⁡(π(t+c)), where a≥0,b and c are constants.By integrating both sides of (16), this paper obtains the exact solution of the instantaneous intensity: (17)ln⁡mt=ln⁡m0+∫0tαsds+∫0tνdWms=ln⁡m0+∫0tαsds+νWms.Then, (18)mt=m0exp⁡∫0tαsds+νWmt,or for all 0≤u≤t(19)mt=muexp⁡∫utαsds+νWmt-Wmu,where αt=a+bsin⁡(π(t+c)),a≥0, andb and c are constants. Figure 6 presents an example of a trajectory followed by the BDT intensity process. It can be observed that these curves exhibit certain periodic and upward trends.Figure 6
Example of path for the stochastic intensitymt. The graph corresponds to the parameter values m0=20, a=0.05, b=-0.3, and c=-0.5.


## 2.4.2. Modeling Loss Severities

Extreme value theory (EVT) provides a simple technique for evaluating probabilities of future extreme levels of a process based on historical observations. Generally, being a large claim means exceeding some prespecified threshold value. The generalized Pareto distribution (GPD) can be used to set the value of this threshold (see for a detailed account in the book by Embrechts et al. [17]). Suppose insured claims are denoted by random variables X over a high threshold u with distribution function F. Then the distribution is defined by(20)Fuy=pX≤u+y∣X>u=Fy+u-Fu1-Fu=Fx-Fu1-Fu,which denotes the probability that the value of x exceeds the threshold u by at most an amount y given that x exceeds the threshold u. Balkema and de Haan [18] and Pickands III [19] showed that, for sufficiently large threshold u, the distribution function of the excess can be approximated by the GPD such that, as the threshold u gets high, the excess distribution Fu(y) converges to the GPD, which is(21)Gξ,σXx=1-1+ξxσX-1/ξifξ≠0,1-exp⁡-xσXifξ=0,where σX>0, and the range for x is x≥0, when ξ≥0, and x∈[0,-σX/ξ] when ξ<0. The parameters σ and ξ are expressed as the scale and shape parameters, respectively.Suppose thatX={X1,X2,…,Xnu} is a random sampling of the GPD, where nu denote the sample exceeding a threshold level u. For each of these exceedances we calculate the amount yi=Xi-u(i=1,2,…,nu) of the excess loss. These generalized Pareto distribution parameters may be evaluated by maximum likelihood estimation (MLE) using the following stages:(1)
Make log-likelihood function forξ≠0(22)lθ∣y=ln⁡Lθ∣y=-nuln⁡σX-∑i=1nu1+1ξln⁡1+ξyiσX,whereθ=(ξ,σX), σX>0, 1+ξyi/σX>0.(2)
Solve the maximization of (22) by using numerical optimization algorithms. Through differentiating (22) with respect to ξ and σX and then equating to zero, this paper gets the following equations: (23)∂lθ∣y∂ξ=1ξ2∑i=1nuln⁡1+ξyiσX+ξyi-∑i=1nuyiσX+ξyi=0,∂lθ∣y∂σX=-nuσX+1ξ+1∑i=1nuξyiσX2+σXξyi=0.


## 2.4.3. Doubly Stochastic Compound Poisson Process

The risk faced by the (re)insurance industry is intrinsic in their exposure to aggregate claims. This paper assumes that the process of aggregate claims is a compound doubly stochastic Poisson process. Then we need to make some assumptions about the aggregate claims process.(i)
There exists a doubly stochastic Poisson processMs(s∈[0,T]) representing potentially catastrophic events. The Poisson point arrival process is assumed to be a predictable bounded process ms. This paper denotes the moments of these potentially catastrophic risk as 0≤t1≤⋯≤⋯≤tj≤⋯≤T.(ii)
The insured claims incurred by each catastrophic risk in the flow{tj}j=1,… are assumed to be independently and identically distributed random variables {Xj}j=1,… with a distribution function F(x)=P{Xj<x} and a density function f(x).Therefore, a left-continuous and predictable aggregate claims process can be defined as(24)Lt=∑tj<tXj=∑j=1MtXj,where Mt and Xj are assumed to be independent.In addition, there is also another assumption:(iii)
The threshold level is the timeτ when the aggregate claims Lt exceeds the threshold value D; then, τ=inf{t:Lt≥D}. Now defining a new process Nt=1{Lt≥D}, Baryshnikov et al. [20] showed that it is a doubly stochastic Poisson process with stochastic intensity: (25)λs=ms1-FD-Ls1Ls<D.However, in that paper, the proof ofλt is rather loose. So, this paper will give a detailed proof of the result in Appendix A.


## 3. Valuation of Zero-Coupon CAT Bonds

In this section, this paper adopts a compound doubly stochastic Poisson process pricing methodology and incorporates more realistic interest rates into models. Before deriving a zero-coupon CAT bond pricing formula, this paper first derives the equivalent martingale measure when the arrival process of catastrophic risk follows a doubly stochastic Poisson process with stochastic intensity.


### 3.1. Equivalent Martingale Probability Measure

Because the natural catastrophe claims process has jumps, the catastrophe bond market is incomplete and there is no unique equivalent martingale measure. Hence, catastrophe bond price could not be unique. This paper adopts Merton’s assumption that the risk associated with jumps can be diversified away [21]. That is, there exits a zero risk premium (this stance is supported by the empirical studies of Cummins and Weiss [22]). Thus, the jump activity rate and distribution are not altered from the physical measure P to the risk-neutral measure Q.Proposition 1.
Letη1(t) denote the Radon-Nikodym derivative process for the CAT process and the Hull-White process as follows: (26)log⁡η1t=log⁡dQdPt=-12∫0tλr2sds+∫0tλrsdWrs+∑j=1Mmtlog⁡φXj,tj+∫0t∫0∞1-φx,smsfxdxds,where ∫0tλr2(u)du<∞ and ∫0t∫0∞[1-φ(x,s)msf(x)]dxdu<∞. φ(x,s),λr(u) may be expressed as the component of the market prices associated with the CAT risk premium associated with the CAT loss and the Winner Process for the interest rate, respectively. For detailed discussion of the jump risk premium φ(x,s), see Glasserman and Kou [23]. By the Radon-Nikodym process, the new Winner Process of the interest rate WrQ(t) is defined by WrQ(t)=Wr(t)+∫0tλr(s)ds. Then the new arrival rate at time s,s∈[0,t] becomes msQ=φ(x,s)msf(x)dx.Note that this paper follows the Merton measure and has a zero risk premium. That is,φ(x,s)=1.The forward martingale measureQT is defined to be the equivalent measure when choosing the T maturity discount bond as the numèraire asset. This paper assumes that the interest rate follows the Hull-White model, which has the following analytical expression:(27)Dt,T=At,Te-Bt,Trt,where (28)Bt,T=1κ1-e-κT-t,log⁡At,T=log⁡D0,TD0,t-Bt,T∂log⁡D0,t∂t-σr24κ1-e-2κtBt,T2.Hence, the price process satisfies (29)dDt,TDt,T=rtdt-σrBt,TdWrQt.Letη2(t) denote the Radon-Nikodym derivative process for forward measure QT as follows: (30)η2t=log⁡dQTdQt=-12∫0tσr2B2s,TdS-∫0tσrBs,TdWsQ,where η2(0)=1 and suppose EQ[η2(t)]=1 for all t. According to Girsanov’s theorem, there exists a forward measure QT such that (31)WrQTt=WrQt+∫0tσrBs,Tds.


### 3.2. The Pricing Formula of Zero-Coupon CAT Bonds

Once the forward-neutral measure is determined, a contingent claim can be priced by the discounted expectation of its various payoffs. This paper considers a zero-coupon CAT bond at maturity timeT contingent on threshold time τ>T. Defining the process PCAT(s)=E(PCAT∣Fs), which requires that PCAT(s) be a predictable process, it can be expressed to mean that the payment at maturity T is not directly linked to the occurrence of the threshold. Therefore, this paper obtains the following result.Theorem 2.
The nonarbitrage price of the zero-coupon CAT bond at timet associated with the threshold level D and the catastrophic flow Ms with stochastic intensity ms, and a loss distribution function F(x) paying principal amount PCAT at time to maturity T under the forward-neutral pricing measure QT is given by (32)VCATt=EQTPCATe-∫tTrsds1-∫tTms1-FD-LS1Ls<Dds∣Ft=BHLt,T,yEPCAT1-∫tTms1-FD-LS1Ls<Dds∣Ft,where (33)BHLt,T,y=exp⁡At,T-Bt,Tyt,At,T=σr22κ3κT-t-21-e-κT-t+121-e-2κT-t-∫tTβsds,Bt,T=1-e-κT-tκ,βT=f0,T+σr22κ21-e-κT2.Proof.
The proof is given in AppendicesB and C.


## 3.1. Equivalent Martingale Probability Measure

Because the natural catastrophe claims process has jumps, the catastrophe bond market is incomplete and there is no unique equivalent martingale measure. Hence, catastrophe bond price could not be unique. This paper adopts Merton’s assumption that the risk associated with jumps can be diversified away [21]. That is, there exits a zero risk premium (this stance is supported by the empirical studies of Cummins and Weiss [22]). Thus, the jump activity rate and distribution are not altered from the physical measure P to the risk-neutral measure Q.Proposition 1.
Letη1(t) denote the Radon-Nikodym derivative process for the CAT process and the Hull-White process as follows: (26)log⁡η1t=log⁡dQdPt=-12∫0tλr2sds+∫0tλrsdWrs+∑j=1Mmtlog⁡φXj,tj+∫0t∫0∞1-φx,smsfxdxds,where ∫0tλr2(u)du<∞ and ∫0t∫0∞[1-φ(x,s)msf(x)]dxdu<∞. φ(x,s),λr(u) may be expressed as the component of the market prices associated with the CAT risk premium associated with the CAT loss and the Winner Process for the interest rate, respectively. For detailed discussion of the jump risk premium φ(x,s), see Glasserman and Kou [23]. By the Radon-Nikodym process, the new Winner Process of the interest rate WrQ(t) is defined by WrQ(t)=Wr(t)+∫0tλr(s)ds. Then the new arrival rate at time s,s∈[0,t] becomes msQ=φ(x,s)msf(x)dx.Note that this paper follows the Merton measure and has a zero risk premium. That is,φ(x,s)=1.The forward martingale measureQT is defined to be the equivalent measure when choosing the T maturity discount bond as the numèraire asset. This paper assumes that the interest rate follows the Hull-White model, which has the following analytical expression:(27)Dt,T=At,Te-Bt,Trt,where (28)Bt,T=1κ1-e-κT-t,log⁡At,T=log⁡D0,TD0,t-Bt,T∂log⁡D0,t∂t-σr24κ1-e-2κtBt,T2.Hence, the price process satisfies (29)dDt,TDt,T=rtdt-σrBt,TdWrQt.Letη2(t) denote the Radon-Nikodym derivative process for forward measure QT as follows: (30)η2t=log⁡dQTdQt=-12∫0tσr2B2s,TdS-∫0tσrBs,TdWsQ,where η2(0)=1 and suppose EQ[η2(t)]=1 for all t. According to Girsanov’s theorem, there exists a forward measure QT such that (31)WrQTt=WrQt+∫0tσrBs,Tds.


## 3.2. The Pricing Formula of Zero-Coupon CAT Bonds

Once the forward-neutral measure is determined, a contingent claim can be priced by the discounted expectation of its various payoffs. This paper considers a zero-coupon CAT bond at maturity timeT contingent on threshold time τ>T. Defining the process PCAT(s)=E(PCAT∣Fs), which requires that PCAT(s) be a predictable process, it can be expressed to mean that the payment at maturity T is not directly linked to the occurrence of the threshold. Therefore, this paper obtains the following result.Theorem 2.
The nonarbitrage price of the zero-coupon CAT bond at timet associated with the threshold level D and the catastrophic flow Ms with stochastic intensity ms, and a loss distribution function F(x) paying principal amount PCAT at time to maturity T under the forward-neutral pricing measure QT is given by (32)VCATt=EQTPCATe-∫tTrsds1-∫tTms1-FD-LS1Ls<Dds∣Ft=BHLt,T,yEPCAT1-∫tTms1-FD-LS1Ls<Dds∣Ft,where (33)BHLt,T,y=exp⁡At,T-Bt,Tyt,At,T=σr22κ3κT-t-21-e-κT-t+121-e-2κT-t-∫tTβsds,Bt,T=1-e-κT-tκ,βT=f0,T+σr22κ21-e-κT2.Proof.
The proof is given in AppendicesB and C.


## 4. Empirical Experiments

This section evaluates zero-coupon CAT bonds whose total losses are linked to the PCS index in the US that occurred between 2001 and 2010 (http://www.verisk.com/insurance/products/property-claim-services/pcs-catastrophe-loss-index.html. A set of 260 catastrophe loss data and 264 catastrophe events in US from 2001 to 2010) and adjust for inflation using the Consumer Price Index (CPI) provided by the US Department of Labor (https://www.bls.gov/cpi/data.htm) That is, the losses are converted to 2010 dollars using the CPI.


### 4.1. Exploratory Data Analysis

There are 260 observations in the data set. Table1 gives a brief summary of descriptive statistics for PCS loss index. According to Table 1, the skewness coefficient is 9.880; then the loss data are considered to be right-skewed. Therefore, this paper utilizes the log-scale to depict the histogram of PCS loss index. In Figure 7 we can discover that the histogram is markedly right-skewed. Table 1 also shows an excess kurtosis of over 115, which is well above the value (i.e., 3). The data are extremely heavy-tailed and should not be modeled by a normal distribution.Table 1
Summary statistics for PCS loss index.
Statistic
Value
Percentile
Value
Sample size
260
Min
0.2541
Range
458.63
5%
0.5471
Mean
8.7048
10%
0.7574
Variance
1213.4
25% (Q1)
1.2697
Std. deviation
34.834
50% (Median)
2.1332
Coef. of variance
4.0017
75% (Q3)
5.0707
Std. error
2.1603
90%
11.451
Skewness
9.8807
95%
24.0
Excess kurtosis
115.24
Max
458.88
Unit: one hundred million dollars.Figure 7
Log histogram of the PCS loss index from 2001 to 2010.To go further to detect the tail behavior, another two graphical techniques deserve special mention: the exponential quantile-quantile (Q-Q) plots and the mean excess function. Figure9 displays the graph of en(u) for PCS loss index. The heavy-tailed character is obvious from the visible upward trend. The concave behavior of the exponential QQ-plot is also apparent in Figure 8, suggesting a heavy-tailed underlying distribution. Therefore, the graphs strongly suggest that the hypotheses that PCS loss index follows a GPD distribution is acceptable.Figure 8
Exponential QQ-plot of the PCS loss index from 2001 to 2010.Figure 9
Mean excess plot of the PCS loss index from 2001 to 2010.


### 4.2. Parameter Estimation of Generalized Pareto Distribution (GPD)

One of the main concerns about the GPD is the selection of an appropriate threshold value. If lower threshold value is applied, the mean exceedance plot could not be linear. To identify the optimal threshold value, two methods are available for this purpose: one is the mean excess function. An immediate result of Section4.1 enables us to provide a graphical tool to choose the threshold value u, through choosing u>0 such that e(x) is approximately linear (tends towards infinity) for x, u. The other is the Hill plot.LetX1>X2>⋯>Xn be the ordered statistics of random observations. The Hill estimator of the tail index ξ using k+1 ordered statistics is defined by (34)Hk,n=1k∑i=1kln⁡XiXk+1,where k is the number of upper order statistics. The Hill plot is the plot of ((k,Hk,n-1),1≤k<n).Based on previous work by Reiss and Thomas [24], Loretan and Phillips [25], and Zhou et al. [26], this paper uses the following rules:(1)
The number of the turning points is not more thank1=[n/10].(2)
The Hill estimator has a relative large deviation.(3)
The turning point is the last number of observation point that satisfies the two above-mentioned conditions.Figure10 shows the Hill plot corresponding to PCS loss index with a 95 percent confidence interval. This paper selects the last area to k=0.1∗260=26, where the Hill plot is more stable. From a closer inspection of the plot in Figure 9, which zooms the function for a smaller range of values for u, this paper suggests selecting the value u=11.49 for the threshold value because it is located at the beginning of a portion of the mean excess plot which is roughly linear. Therefore, this chooses u=11.49.Figure 10
Hill plot versusk upper order statistics (k≤260).Applying the maximum likelihood estimation method, this paper evaluates the valuesξ^=0.9238 and σ^=15.5797 which maximize the log-likelihood function l=-116.7433 for the PCS loss index. Again it may be best to apply a QQ-plot to visually check whether the observations meet the GPD assumption. Figure 11 depicts the plot of the PCS loss index quantiles against Gξ^,σ^ quantiles. So it can be concluded that the fit is satisfactory.Figure 11
QQ-plot of PCS loss index quantiles againstGξ^,σ^ quantiles.


### 4.3. Parameter Estimation of Stochastic Intensity

The parameters of BDT stochastic intensity for the arrival rate of catastrophic events will be estimated. This paper takes the same approach as Lin et al. [27]; given ν=0, by applying the Nonlinear Least Square procedure we conclude that m0=22.3512,a=0.0457,b=0.7338,c=-0.6423; that is, mt=22.3512exp⁡[∫0t0.0457+0.7338sin⁡π(s-0.6423)ds]. Lin et al. [27] proposed a lognormal intensity function of the form λd1(t)=λd10exp⁡[αt-1/2ν2+νWm(t)]. In the same way, it can be concluded that λd1(0)=19.6681,α=0.0515; that is, λd1(t)=19.6681exp⁡[0.0515t]. This paper also considers a pure Poisson process, Burnecki et al. [28] proposed a nonhomogeneous Poisson intensity function of the form λd2(t)=a+b·2π·sin⁡{2π(t-c)}. Using the Nonlinear Least Square procedure we conclude that a=20.1985,b=1.0397 and c=0.1991.


#### 4.3.1. Model Evaluation

To assess these Poisson intensities precisely, it is necessary to calculate the errors based on the comparison of observation value with the values calculated by these Poisson intensities. This paper uses six goodness-of-fit measures to evaluate the performance of the intensity: the root mean square error, the mean absolute error, the mean relative percentage error, Theil’s coefficient, Nash-Sutcliffe model coefficient of efficiency, and the index of agreement:(35)RMSE=1N∑i=1NOi-Pi2,MAE=1N∑i=1NOi-Pi,MRPE=1N∑i=1NOi-PiOi,U=1/N∑i=1NOi-Pi21/N∑i=1NOi2+1/N∑i=1NPi2,E=1.0-∑i=1NOi-Pi2∑i=1NOi-O-2,d=1.0-∑i=1NOi-Pi2∑i=1NPi-O-+Oi-O-2,where Oi denotes the observed value, the overbar represents the mean for the full time period of the evaluation, Pi is the theoretical value, and N is all the number of observations. The MAE, RMSE, and MRPE statistics have, as the lower limit, the value of zero, which is the optimum value for them as it for U. On the contrary, higher values of E and d represent better agreement between the theoretical values and observations [29].Table2 gives a summary of six measurements. Obviously, mt is superior to two other intensities ascribable to its smaller RMSE, MAE, MRPE, and U and larger E and d.Table 2
Summary of the statistical tests for the Poisson intensity.
Poisson intensity
Performance
RMSE
MAE
MAPE
U
E
d
mt
2.3807
1.8595
0.0695
0.0443
0.7964
0.9406
λ
d
1
(
t
) 3.5495
2.9540
0.1099
0.0602
0.5475
0.8341
λd2(t)
5.2764
4.6800
0.1780
0.0990
-2.5791e-011
9.1819e-006
Note: the optimum value is in boldface.


## 4.1. Exploratory Data Analysis

There are 260 observations in the data set. Table1 gives a brief summary of descriptive statistics for PCS loss index. According to Table 1, the skewness coefficient is 9.880; then the loss data are considered to be right-skewed. Therefore, this paper utilizes the log-scale to depict the histogram of PCS loss index. In Figure 7 we can discover that the histogram is markedly right-skewed. Table 1 also shows an excess kurtosis of over 115, which is well above the value (i.e., 3). The data are extremely heavy-tailed and should not be modeled by a normal distribution.Table 1
Summary statistics for PCS loss index.
Statistic
Value
Percentile
Value
Sample size
260
Min
0.2541
Range
458.63
5%
0.5471
Mean
8.7048
10%
0.7574
Variance
1213.4
25% (Q1)
1.2697
Std. deviation
34.834
50% (Median)
2.1332
Coef. of variance
4.0017
75% (Q3)
5.0707
Std. error
2.1603
90%
11.451
Skewness
9.8807
95%
24.0
Excess kurtosis
115.24
Max
458.88
Unit: one hundred million dollars.Figure 7
Log histogram of the PCS loss index from 2001 to 2010.To go further to detect the tail behavior, another two graphical techniques deserve special mention: the exponential quantile-quantile (Q-Q) plots and the mean excess function. Figure9 displays the graph of en(u) for PCS loss index. The heavy-tailed character is obvious from the visible upward trend. The concave behavior of the exponential QQ-plot is also apparent in Figure 8, suggesting a heavy-tailed underlying distribution. Therefore, the graphs strongly suggest that the hypotheses that PCS loss index follows a GPD distribution is acceptable.Figure 8
Exponential QQ-plot of the PCS loss index from 2001 to 2010.Figure 9
Mean excess plot of the PCS loss index from 2001 to 2010.


## 4.2. Parameter Estimation of Generalized Pareto Distribution (GPD)

One of the main concerns about the GPD is the selection of an appropriate threshold value. If lower threshold value is applied, the mean exceedance plot could not be linear. To identify the optimal threshold value, two methods are available for this purpose: one is the mean excess function. An immediate result of Section4.1 enables us to provide a graphical tool to choose the threshold value u, through choosing u>0 such that e(x) is approximately linear (tends towards infinity) for x, u. The other is the Hill plot.LetX1>X2>⋯>Xn be the ordered statistics of random observations. The Hill estimator of the tail index ξ using k+1 ordered statistics is defined by (34)Hk,n=1k∑i=1kln⁡XiXk+1,where k is the number of upper order statistics. The Hill plot is the plot of ((k,Hk,n-1),1≤k<n).Based on previous work by Reiss and Thomas [24], Loretan and Phillips [25], and Zhou et al. [26], this paper uses the following rules:(1)
The number of the turning points is not more thank1=[n/10].(2)
The Hill estimator has a relative large deviation.(3)
The turning point is the last number of observation point that satisfies the two above-mentioned conditions.Figure10 shows the Hill plot corresponding to PCS loss index with a 95 percent confidence interval. This paper selects the last area to k=0.1∗260=26, where the Hill plot is more stable. From a closer inspection of the plot in Figure 9, which zooms the function for a smaller range of values for u, this paper suggests selecting the value u=11.49 for the threshold value because it is located at the beginning of a portion of the mean excess plot which is roughly linear. Therefore, this chooses u=11.49.Figure 10
Hill plot versusk upper order statistics (k≤260).Applying the maximum likelihood estimation method, this paper evaluates the valuesξ^=0.9238 and σ^=15.5797 which maximize the log-likelihood function l=-116.7433 for the PCS loss index. Again it may be best to apply a QQ-plot to visually check whether the observations meet the GPD assumption. Figure 11 depicts the plot of the PCS loss index quantiles against Gξ^,σ^ quantiles. So it can be concluded that the fit is satisfactory.Figure 11
QQ-plot of PCS loss index quantiles againstGξ^,σ^ quantiles.


## 4.3. Parameter Estimation of Stochastic Intensity

The parameters of BDT stochastic intensity for the arrival rate of catastrophic events will be estimated. This paper takes the same approach as Lin et al. [27]; given ν=0, by applying the Nonlinear Least Square procedure we conclude that m0=22.3512,a=0.0457,b=0.7338,c=-0.6423; that is, mt=22.3512exp⁡[∫0t0.0457+0.7338sin⁡π(s-0.6423)ds]. Lin et al. [27] proposed a lognormal intensity function of the form λd1(t)=λd10exp⁡[αt-1/2ν2+νWm(t)]. In the same way, it can be concluded that λd1(0)=19.6681,α=0.0515; that is, λd1(t)=19.6681exp⁡[0.0515t]. This paper also considers a pure Poisson process, Burnecki et al. [28] proposed a nonhomogeneous Poisson intensity function of the form λd2(t)=a+b·2π·sin⁡{2π(t-c)}. Using the Nonlinear Least Square procedure we conclude that a=20.1985,b=1.0397 and c=0.1991.


### 4.3.1. Model Evaluation

To assess these Poisson intensities precisely, it is necessary to calculate the errors based on the comparison of observation value with the values calculated by these Poisson intensities. This paper uses six goodness-of-fit measures to evaluate the performance of the intensity: the root mean square error, the mean absolute error, the mean relative percentage error, Theil’s coefficient, Nash-Sutcliffe model coefficient of efficiency, and the index of agreement:(35)RMSE=1N∑i=1NOi-Pi2,MAE=1N∑i=1NOi-Pi,MRPE=1N∑i=1NOi-PiOi,U=1/N∑i=1NOi-Pi21/N∑i=1NOi2+1/N∑i=1NPi2,E=1.0-∑i=1NOi-Pi2∑i=1NOi-O-2,d=1.0-∑i=1NOi-Pi2∑i=1NPi-O-+Oi-O-2,where Oi denotes the observed value, the overbar represents the mean for the full time period of the evaluation, Pi is the theoretical value, and N is all the number of observations. The MAE, RMSE, and MRPE statistics have, as the lower limit, the value of zero, which is the optimum value for them as it for U. On the contrary, higher values of E and d represent better agreement between the theoretical values and observations [29].Table2 gives a summary of six measurements. Obviously, mt is superior to two other intensities ascribable to its smaller RMSE, MAE, MRPE, and U and larger E and d.Table 2
Summary of the statistical tests for the Poisson intensity.
Poisson intensity
Performance
RMSE
MAE
MAPE
U
E
d
mt
2.3807
1.8595
0.0695
0.0443
0.7964
0.9406
λ
d
1
(
t
) 3.5495
2.9540
0.1099
0.0602
0.5475
0.8341
λd2(t)
5.2764
4.6800
0.1780
0.0990
-2.5791e-011
9.1819e-006
Note: the optimum value is in boldface.


## 4.3.1. Model Evaluation

To assess these Poisson intensities precisely, it is necessary to calculate the errors based on the comparison of observation value with the values calculated by these Poisson intensities. This paper uses six goodness-of-fit measures to evaluate the performance of the intensity: the root mean square error, the mean absolute error, the mean relative percentage error, Theil’s coefficient, Nash-Sutcliffe model coefficient of efficiency, and the index of agreement:(35)RMSE=1N∑i=1NOi-Pi2,MAE=1N∑i=1NOi-Pi,MRPE=1N∑i=1NOi-PiOi,U=1/N∑i=1NOi-Pi21/N∑i=1NOi2+1/N∑i=1NPi2,E=1.0-∑i=1NOi-Pi2∑i=1NOi-O-2,d=1.0-∑i=1NOi-Pi2∑i=1NPi-O-+Oi-O-2,where Oi denotes the observed value, the overbar represents the mean for the full time period of the evaluation, Pi is the theoretical value, and N is all the number of observations. The MAE, RMSE, and MRPE statistics have, as the lower limit, the value of zero, which is the optimum value for them as it for U. On the contrary, higher values of E and d represent better agreement between the theoretical values and observations [29].Table2 gives a summary of six measurements. Obviously, mt is superior to two other intensities ascribable to its smaller RMSE, MAE, MRPE, and U and larger E and d.Table 2
Summary of the statistical tests for the Poisson intensity.
Poisson intensity
Performance
RMSE
MAE
MAPE
U
E
d
mt
2.3807
1.8595
0.0695
0.0443
0.7964
0.9406
λ
d
1
(
t
) 3.5495
2.9540
0.1099
0.0602
0.5475
0.8341
λd2(t)
5.2764
4.6800
0.1780
0.0990
-2.5791e-011
9.1819e-006
Note: the optimum value is in boldface.


## 5. Numerical Analysis

This section considers a numerical analysis for the prices of zero-coupon CAT bonds using the Monte Carlo method with 20000 paths. First, the parameter of the stochastic interest rates applied in our numerical analysis were fitted in Hull and White [10]. Then, a base set of parameters is summarized in Table 3.Table 3
Parameters, definitions, and base values.
Parameters
Description
Values
Interest rate parameters
κ
Magnitude of mean-reverting force
0.1
σ
r
Volatility of interest rate
0.01, 0.03, 0.05
Z
t
Continuously compounded zero-coupon rate
0.08-0.05e-0.18t
Catastrophe loss parameters
ξX
Shape parameter of the GPD of the amount of catastrophe losses for the insurer
0.9238
σ
X
Scale parameter of the GPD of the amount of catastrophe losses for the insurer
[12,16]
a
,
b
,
c
Drift parameters of intensity function for catastrophic events
0.0457,0.7388,-0.6423
ν
Volatility of intensity function for catastrophic events
[
0.5,1.0
]
m
0
Initial value of intensity function for catastrophic events
22.3512
Other parameters
T
Time to expiry
2
t
Trading time
[
0,2
]
V
Face value of the CAT bond
1
D
Threshold levels
[10000,20000]Figure12 shows that the zero-coupon CAT bond prices with respect to trading time t to expiry T and threshold level D under the GPD, DSPP, and Hull-White interest rate assumptions. It is obvious that the zero-coupon CAT bond prices increase with the trading time t as the occurrence probability exceeding a threshold value will decrease, and it decreases as the threshold level D decreases because one expects a high probability exceeding threshold value. This means that the increase in trading time t reduces the probability of a trigger event, causing a higher price, and simultaneously much less coupons are expected to be received, which also leads to a higher price. When t=1.95 years and D=20000, the zero-coupon CAT bond price approaches the value exp⁡{A(1.95,2)-B(1.95,2)y1.9}≈0.99659.Figure 12
The zero-coupon CAT bond price with respect to the threshold level and trading time. The model parameters areξX=0.9238,σX=15.5797,m0=22.3512,a=0.0457,b=0.7388,c=-0.6423,ν=0.5,κ=0.1, and σr=0.01.To measure how stochastic interest rates influence the zero-coupon CAT bond prices, this paper considers the price difference from the Hull-White model and a discount bond. When timet=0, the Hull-White extended-Vasicek model value is equal with the discount bond D(0,T)=exp{-ZTT}, when time t=T, BHL(T,T,yT)=D(T,T)=1; that is, D(t,T)=exp{-ZT-t×(T-t)}. Figure 13 shows the price difference between Hull-White interest rates and the deterministic interest rate. It can be observed that the price differences increase first and subsequently decrease, and the deterministic interest rate overestimates the bond price change substantially from trading time 0 to 2 years.Figure 13
The difference between the zero-coupon CAT bond price in stochastic interest rates and the deterministic interest rate under the lognormal and the DSPP assumptions. The model parameters areξX=0.9238,σX=15.5797,m0=22.3512,a=0.0457,b=0.7388,c=-0.6423,ν=0.5,κ=0.1, and σr=0.01.Finally, this paper presents a sensitivity analysis for the zero-coupon CAT bond prices with different parameter values using the DSPP and the GPD. Table4 displays the zero-coupon CAT bond prices under various parameter values of the threshold level, the CAT intensity volatility, and loss scale parameter. As is expected, the higher the threshold level, the smaller the CAT intensity volatility and CAT loss scale parameter and the higher the zero-coupon CAT bond prices. For example, when the CAT intensity volatility (ν) is 0.8 and the loss scale parameter (σX) is 13 and the zero-coupon CAT bond price will increase by 519 basis points as the threshold value increases from 14000 to 15000, while the zero-coupon CAT bond price will increase by 416 basis points as the threshold value increases from 15000 to 16000.Table 4
The zero-coupon CAT bond prices with trading timet to expiry T.
(
σ
X
,
ν
)
Threshold level (D)
10000
12000
14000
15000
16000
18000
200000
(11,0.5)
0.85979
0.87052
0.87703
0.87940
0.88139
0.88451
0.88686
(11,0.6)
0.85291
0.86607
0.87370
0.87644
0.87871
0.88225
0.88488
(11,0.7)
0.84807
0.86383
0.87277
0.87595
0.87857
0.88263
0.88563
(11,0.8)
0.82802
0.84938
0.86050
0.86432
0.86744
0.87443
0.87900
(11,0.9)
0.79958
0.83459
0.85176
0.85713
0.86136
0.86765
0.87211
(11,1.0)
0.71016
0.78851
0.83038
0.84227
0.85026
0.86082
0.86766
(12,0.5)
0.85435
0.86670
0.87409
0.87677
0.87907
0.88252
0.88514
(12,0.6)
0.84583
0.86149
0.87028
0.87341
0.87599
0.87999
0.88296
(12,0.7)
0.83958
0.85842
0.86888
0.87245
0.87544
0.88005
0.88345
(12,0.8)
0.81497
0.84244
0.85567
0.86014
0.86374
0.86921
0.87583
(12,0.9)
0.77814
0.82329
0.84477
0.85132
0.85638
0.86376
0.86892
(12,1.0)
0.68427
0.76604
0.81381
0.82973
0.84100
0.85456
0.86289
(13,0.5)
0.84853
0.86269
0.87103
0.87405
0.87655
0.88046
0.88338
(13,0.6)
0.83771
0.85659
0.86670
0.87025
0.87316
0.87766
0.88098
(13,0.7)
0.82996
0.85261
0.86462
0.86878
0.87218
0.87738
0.88119
(13,0.8)
0.80197
0.83468
0.85048
0.85567
0.85983
0.86607
0.87155
(13,0.9)
0.75495
0.81135
0.83677
0.84487
0.85094
0.85960
0.86555
(13,1.0)
0.66265
0.73558
0.79589
0.81432
0.82917
0.84737
0.85762
(14,0.5)
0.84230
0.85849
0.86787
0.87123
0.87401
0.87835
0.88158
(14,0.6)
0.82974
0.85130
0.86293
0.86694
0.87022
0.87525
0.87894
(14,0.7)
0.81878
0.84631
0.86019
0.86493
0.86876
0.87461
0.87886
(14,0.8)
0.78749
0.82580
0.84486
0.85090
0.85567
0.86278
0.86783
(14,0.9)
0.73149
0.79677
0.82715
0.83756
0.84494
0.85514
0.86199
(14,1.0)
0.64217
0.70692
0.77544
0.79836
0.81474
0.83888
0.85173
(15,0.5)
0.83556
0.85407
0.86458
0.86831
0.87140
0.87619
0.87973
(15,0.6)
0.82113
0.84553
0.85894
0.86348
0.86715
0.87275
0.87684
(15,0.7)
0.80613
0.83943
0.85549
0.86086
0.86518
0.87173
0.87645
(15,0.8)
0.77167
0.81523
0.83871
0.84576
0.85125
0.85932
0.86499
(15,0.9)
0.70721
0.77951
0.81801
0.82894
0.83823
0.85032
0.85820
(15,1.0)
0.62390
0.68638
0.75536
0.77876
0.80034
0.82825
0.84502
(16,0.5)
0.82819
0.84942
0.86115
0.86530
0.86870
0.87396
0.87784
(16,0.6)
0.81226
0.83908
0.85472
0.85983
0.86394
0.87017
0.87467
(16,0.7)
0.79374
0.83183
0.85047
0.85657
0.86143
0.86873
0.87396
(16,0.8)
0.75732
0.80475
0.83192
0.84020
0.84652
0.85567
0.86202
(16,0.9)
0.68834
0.76069
0.80744
0.82038
0.83045
0.84507
0.85416
(16,1.0)
0.60683
0.66721
0.72777
0.76147
0.78190
0.81540
0.83718
Note: other parameter values areξX=0.9238,m0=22.3512,a=0.0457,b=0.7388,c=-0.6423,κ=0.1 and σr=0.05, t=0.8, and T=2.


## 6. Conclusion

Consider an upward trend-cycle movement of catastrophic risk events and the severity of aggregate catastrophe claims; this paper reinvestigates and incorporates these catastrophic characteristics into the valuation of zero-coupon CAT bonds and provides a semianalytical solution for evaluating zero-coupon catastrophe bonds with the compound doubly stochastic Poisson process pricing framework under the forward measure. Using the PCS loss data from 2001 to 2010, the empirical results reveal that the PCS loss index following a GPD distribution is acceptable and the fit is satisfactory, and the BDT arrival rate is more superior to the mean-reverting arrival process and the deterministic time intensity due to its smaller RMSE, MAE, MRPE, and U and larger E and d. This paper briefly demonstrates how catastrophic components, the CAT intensity volatility, loss scale parameter, and financial components, as well as the interest rate volatility and the threshold level affect the value of zero-coupon CAT bonds. In short, the numerical results show that the proposed model is feasible and effective.



---

*Source: 3279647-2017-09-13.xml*