Research Article  Open Access
Pricing ZeroCoupon Catastrophe Bonds Using EVT with Doubly Stochastic Poisson Arrivals
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 quantilequantile (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 zerocoupon 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 zerocoupon catastrophe bonds.
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 Figures 1, 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].
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 zerocoupon CAT bonds with a doubly stochastic compound Poisson process pricing framework under the forward measure.
The paper is organized as follows. Section 2 describes the model assumptions and illustrates CAT claim models. In Section 3 we show the forwardneutral measure and derive the pricing formula for zerocoupon 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 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 denote the aggregate loss process of the insured up to time . is a doubly stochastic Poisson point process with an intensity parameter . is a sequence of identically and independently distributed (i.i.d.) random variables and accounts for the sizes of CAT loss.
In addition, let stand for a doubly stochastic Poisson process with stochastic intensity and represent the CAT number, and let denote the riskfree short rate process (or the force of interest). and 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 define a probability space, where is the set of states of the world, the filtration is algebra of subsets of and is defined by , where is and is a probability measure on .
Let denote the price at time of a (defaultfree) discount bond that pays USD 1 at maturity date , which is given byThen, the spot forward rate at time for future date is defined by
2.2. Valuation Theory
In a riskfree arbitrage financial market, the price of the contingent claim at time based on an underlying asset at maturity is given bywhere denotes the conditional expectation and 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 forwardneutral measure . 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 maturity bond as a numeraire for contingent claims that have a payoff at time . By applying the Change of Numeraire Theorem, the prices are martingales for under the measure . Hence, by applying the definition of a martingale, we obtainHowever, at maturity , the price of the discount bond and the price of the contingent claim are given by the payoff . So, (4) reduces toTaking derivatives of (1) with respect to yieldswhere we have used the Change of Numeraire Theorem in the last step with RandonNikodym derivativeThen we can simplify (2) to
Therefore, the spot forward rate is a martingale under the measure .
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 singlefactor 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 HullWhite model [10]. The spot interest rate process is denoted by the following stochastic differential equation (SDE): where and are positive constants and is a deterministic function of time. The parameter here can be analytically computed as
Let and . Then, the mean and variance of are given by for .
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 time is a (non)homogeneous Poisson process, denoted by , with the intensity of the counting process varying with time (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 be a point process adapted to , and let be a nonnegative process. Suppose that is measurable, , and that If for all and then is called doubly stochastic Poisson process with intensity .
Equation (13) yields, for all and , where denotes a factorial.
Thus, the probability of no occurrences within the interval for the process with intensity is given by
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: where is a deterministic function of and is a positive constant.
In addition, to describe a complex periodic characteristic of catastrophic risk events, this paper takes here , where and are constants.
By integrating both sides of (16), this paper obtains the exact solution of the instantaneous intensity: Then, or for all where , and and 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.
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 over a high threshold with distribution function . Then the distribution is defined bywhich denotes the probability that the value of exceeds the threshold by at most an amount given that exceeds the threshold . Balkema and de Haan [18] and Pickands III [19] showed that, for sufficiently large threshold , the distribution function of the excess can be approximated by the GPD such that, as the threshold gets high, the excess distribution converges to the GPD, which iswhere , and the range for is , when , and when . The parameters and are expressed as the scale and shape parameters, respectively.
Suppose that is a random sampling of the GPD, where denote the sample exceeding a threshold level . For each of these exceedances we calculate the amount of the excess loss. These generalized Pareto distribution parameters may be evaluated by maximum likelihood estimation (MLE) using the following stages:(1)Make loglikelihood function for â€‰where , , .(2)Solve the maximization of (22) by using numerical optimization algorithms. Through differentiating (22) with respect to and and then equating to zero, this paper gets the following equations:
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 process representing potentially catastrophic events. The Poisson point arrival process is assumed to be a predictable bounded process . This paper denotes the moments of these potentially catastrophic risk as (ii)The insured claims incurred by each catastrophic risk in the flow are assumed to be independently and identically distributed random variables with a distribution function and a density function .
Therefore, a leftcontinuous and predictable aggregate claims process can be defined as where and are assumed to be independent.
In addition, there is also another assumption:(iii)The threshold level is the time when the aggregate claims exceeds the threshold value ; then, Now defining a new process , Baryshnikov et al. [20] showed that it is a doubly stochastic Poisson process with stochastic intensity:
However, in that paper, the proof of is rather loose. So, this paper will give a detailed proof of the result in Appendix A.
3. Valuation of ZeroCoupon 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 zerocoupon 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 to the riskneutral measure .
Proposition 1. Let denote the RadonNikodym derivative process for the CAT process and the HullWhite process as follows: where and . 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 , see Glasserman and Kou [23]. By the RadonNikodym process, the new Winner Process of the interest rate is defined by . Then the new arrival rate at time becomes .
Note that this paper follows the Merton measure and has a zero risk premium. That is, .
The forward martingale measure is defined to be the equivalent measure when choosing the maturity discount bond as the numÃ¨raire asset. This paper assumes that the interest rate follows the HullWhite model, which has the following analytical expression:where Hence, the price process satisfies
Let denote the RadonNikodym derivative process for forward measure as follows: where and suppose for all . According to Girsanovâ€™s theorem, there exists a forward measure such that
3.2. The Pricing Formula of ZeroCoupon CAT Bonds
Once the forwardneutral measure is determined, a contingent claim can be priced by the discounted expectation of its various payoffs. This paper considers a zerocoupon CAT bond at maturity time contingent on threshold time . Defining the process , which requires that be a predictable process, it can be expressed to mean that the payment at maturity is not directly linked to the occurrence of the threshold. Therefore, this paper obtains the following result.
Theorem 2. The nonarbitrage price of the zerocoupon CAT bond at time associated with the threshold level and the catastrophic flow with stochastic intensity , and a loss distribution function paying principal amount at time to maturity under the forwardneutral pricing measure is given by where
Proof. The proof is given in Appendices B and C.
4. Empirical Experiments
This section evaluates zerocoupon 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/propertyclaimservices/pcscatastrophelossindex.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. Table 1 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 rightskewed. Therefore, this paper utilizes the logscale to depict the histogram of PCS loss index. In Figure 7 we can discover that the histogram is markedly rightskewed. Table 1 also shows an excess kurtosis of over 115, which is well above the value (i.e., 3). The data are extremely heavytailed and should not be modeled by a normal distribution.
 
Unit: one hundred million dollars. 
To go further to detect the tail behavior, another two graphical techniques deserve special mention: the exponential quantilequantile (QQ) plots and the mean excess function. Figure 9 displays the graph of for PCS loss index. The heavytailed character is obvious from the visible upward trend. The concave behavior of the exponential QQplot is also apparent in Figure 8, suggesting a heavytailed underlying distribution. Therefore, the graphs strongly suggest that the hypotheses that PCS loss index follows a GPD distribution is acceptable.
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 Section 4.1 enables us to provide a graphical tool to choose the threshold value , through choosing such that is approximately linear (tends towards infinity) for , . The other is the Hill plot.
Let be the ordered statistics of random observations. The Hill estimator of the tail index using ordered statistics is defined by where is the number of upper order statistics. The Hill plot is the plot of .
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 than .(2)The Hill estimator has a relative large deviation.(3)The turning point is the last number of observation point that satisfies the two abovementioned conditions.
Figure 10 shows the Hill plot corresponding to PCS loss index with a 95 percent confidence interval. This paper selects the last area to , 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 , this paper suggests selecting the value 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 .
Applying the maximum likelihood estimation method, this paper evaluates the values and which maximize the loglikelihood function for the PCS loss index. Again it may be best to apply a QQplot to visually check whether the observations meet the GPD assumption. Figure 11 depicts the plot of the PCS loss index quantiles against quantiles. So it can be concluded that the fit is satisfactory.
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 , by applying the Nonlinear Least Square procedure we conclude that ; that is, . Lin et al. [27] proposed a lognormal intensity function of the form . In the same way, it can be concluded that ; that is, . This paper also considers a pure Poisson process, Burnecki et al. [28] proposed a nonhomogeneous Poisson intensity function of the form . Using the Nonlinear Least Square procedure we conclude that and .
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 goodnessoffit 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, NashSutcliffe model coefficient of efficiency, and the index of agreement:where denotes the observed value, the overbar represents the mean for the full time period of the evaluation, is the theoretical value, and 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 . On the contrary, higher values of E and d represent better agreement between the theoretical values and observations [29].
Table 2 gives a summary of six measurements. Obviously, is superior to two other intensities ascribable to its smaller RMSE, MAE, MRPE, and U and larger E and d.
 
Note: the optimum value is in boldface. 
5. Numerical Analysis
This section considers a numerical analysis for the prices of zerocoupon 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.

Figure 12 shows that the zerocoupon CAT bond prices with respect to trading time to expiry and threshold level under the GPD, DSPP, and HullWhite interest rate assumptions. It is obvious that the zerocoupon CAT bond prices increase with the trading time as the occurrence probability exceeding a threshold value will decrease, and it decreases as the threshold level decreases because one expects a high probability exceeding threshold value. This means that the increase in trading time 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 years and , the zerocoupon CAT bond price approaches the value .
To measure how stochastic interest rates influence the zerocoupon CAT bond prices, this paper considers the price difference from the HullWhite model and a discount bond. When time , the HullWhite extendedVasicek model value is equal with the discount bond , when time , ; that is, . Figure 13 shows the price difference between HullWhite 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.
Finally, this paper presents a sensitivity analysis for the zerocoupon CAT bond prices with different parameter values using the DSPP and the GPD. Table 4 displays the zerocoupon 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 zerocoupon CAT bond prices. For example, when the CAT intensity volatility is 0.8 and the loss scale parameter is 13 and the zerocoupon CAT bond price will increase by 519 basis points as the threshold value increases from 14000 to 15000, while the zerocoupon CAT bond price will increase by 416 basis points as the threshold value increases from 15000 to 16000.
 
Note: other parameter values are and , , and . 
6. Conclusion
Consider an upward trendcycle movement of catastrophic risk events and the severity of aggregate catastrophe claims; this paper reinvestigates and incorporates these catastrophic characteristics into the valuation of zerocoupon CAT bonds and provides a semianalytical solution for evaluating zerocoupon 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 meanreverting 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 zerocoupon CAT bonds. In short, the numerical results show that the proposed model is feasible and effective.
Appendix
A. Proof of (25)
This appendix is devoted to the proof of doubly stochastic Poisson process with intensity .
Representing as and from the definition of , it follows thatwhere denotes the index of the first event following (finite a.e.); that is, .
Let be a nonnegative predictable process. Then, by using the basic properties of conditional expectation and Fubiniâ€™s theorem, we haveFor the last equality, use the fact that is measurable. Since is independent of , we haveLet . Therefore, we have which is verified for all nonnegative predictable process ; then is a doubly stochastic Poisson process with intensity .
B. The Price of a Discount Bond
The price of the contingent claim at time , which has a payoff at maturity satisfies the partial differential equation: where follows the HullWhite model and denotes the market price of interest rate risk.
First, consider the following transformation of variables. The process can be written as follows: where , which guarantees .
The price of the contingent claim in terms of the new variable can be written as ; that is, . We can establish the following relationships between and :Substituting these relations into (B.1) and using , the partial differential equation reduces to Under the equivalent martingale measure , the spot interest rate is generated bywhere is still a Wiener process.
Given a value at any point in time , the probability distribution of for some future time is a normal distribution with mean, and variance,
Using the FeynmanKac formula, we can express solutions to the partial differential equation (B.3) in terms of the boundary condition at time as Applying the forward measure , the price can be expressed as where denotes the price of a discount bond with maturity at time and is given by