Pricing Multi-Event-Triggered Catastrophe Bonds Based on a Copula–POT Model

Abstract

The constantly expanding losses caused by frequent natural disasters pose many challenges to the traditional catastrophe insurance market. The purpose of this paper is to develop an innovative and systemic trigger mechanism for pricing catastrophic bonds triggered by multiple events with an extreme dependence structure. Due to the bond’s low cashflow contingencies and the CAT bond’s high return, the multiple-event CAT bond may successfully transfer the catastrophe risk to the huge financial markets to meet the diversification of capital allocations for most potential investors. The designed hybrid trigger mechanism helps reduce the moral hazard and increase the bond’s attractiveness with a lower trigger likelihood, displaying the determinants of the wiped-off coupon and principal by both the magnitude and intensity of the natural disaster events involved. As the trigger indicators resulting from the potential catastrophic disaster might be associated with heavy-tailed margins, nested Archimedean copulas are introduced with marginal distributions modeled by a POT-GP distribution for excess data and common parametric models for moderate risks. To illustrate our theoretical pricing framework, we conduct an empirical analysis of pricing a three-event rainstorm CAT bond based on the resulting losses due to rainstorms in China during 2006–2020. Monte Carlo simulations are carried out to analyze the sensitivity of the rainstorm CAT bond price in trigger attachment levels, maturity date, catastrophe intensity, and numbers of trigger indicators.

1. Introduction

In recent years, the deterioration of the natural environment and growing human activities have increased the level of damage caused by natural disasters, and the economic losses incurred have been on the rise. According to the Sigma Catastrophe Database, the global catastrophe data during 1970–2021 show that the frequency of catastrophes has generally risen over the past fifty years. The same goes for losses, with 90% of all losses over the last decade being in the tens of billions of dollars or more, placing a heavy burden on insurance companies, government finances, and society.

Many countries disperse the catastrophe risk by issuing various types of catastrophe securitization products. For example, in 1997, Hannover Re launched the first successful issue of a CAT bond that included exposure to hurricane and earthquake disasters in Japan, Australia, and Canada. This form filled in a gap in the traditional insurance approach, transferring disaster risks to the huge capital market, effectively diversifying catastrophe risk and improving the payment capacity of insurance companies (Braun 2016; Karagiannis et al. 2016). Among these, catastrophe bonds are considered to be the most mature financial instrument among catastrophe securitization products, attracting an increasing number of investors because it is almost uncorrelated with the returns of other financial market instruments (Richter and Wilson 2020).

A CAT bond is a financial security that pays the issuer (i.e., collateralized special purpose vehicles (SPVs)) when a predefined disaster risk materializes (Tang and Yuan 2019). Figure 1 illustrates its flowchart. These SPVs receive reinsurance premiums from investors and re-insurers in exchange for providing reinsurance coverage. The premiums, often part of coupon payments, may also include a floating portion linked to a reference rate like LIBOR, reflecting the trustee’s return. Upon the occurrence of the specified triggering event, the principal and coupon payments are reduced to allocate funds for reimbursing the sponsors for claims paid (Cummins 2008). Notably, the trigger mechanism, determining when payouts are made to the bond issuer, is a key distinguishing feature among issuer types (Cox and Pedersen 2000).

As CAT bonds are designed to cover high layers of insurance losses, a single unprecedented event caused by catastrophic disasters might be selected as a trigger event (Aase 1999; Chen et al. 2013; Litzenberger et al. 1996; Nowak and Romaniuk 2013). Although such industry loss trigger mechanisms can completely eliminate the basis risk for insurers, the moral hazard is very high, since insurers may exaggerate losses in their loss statistics for their own benefit, generating a vast information asymmetry with capital market investors (Pauly 1968). Additionally, as argued by Ibrahim et al. (2022), investors’ interest in single-event-triggered catastrophe bonds is likely to decline in the future, since the global trend of increasing year-on-year disaster losses and intensity may increase the risk of catastrophe bond claims. Thus, issuing multi-trigger catastrophe bonds could be a solution using both industry loss indices and physical parameters as the triggering conditions (Cox and Pedersen 2000; Woo 2004). This hybrid trigger mechanism can avoid a basis risk and moral hazard with a reduced triggering risk, attracting investors with a low risk appetite in the market. In addition to understanding the tail behavior of the marginal trigger indicators within the framework of extreme value theory (Anantapadmanabhan 1971; Deng et al. 2020; Leppisaari 2016; McNeil 1997; Zimbidis et al. 2007), the dependence among multiple trigger indicators needs to be considered through a copula approach (Chao and Zou 2018; Reshetar 2008; Wei et al. 2022).

This paper aims to develop an innovative and systemic trigger mechanism in the design of CAT bond pricing models. It allows for a step-wise coupon paid out based on both the intensity of catastrophic events and the severity of its resulting multiple dependent catastrophe indicators. The innovation of our trigger mechanism is four-fold.

• First, the considerable selection of catastrophe indicators, incorporating the hazard information from multiple recording indicators, may result in smaller trigger errors with integrated advantages of the hybrid trigger mechanism. The dependence structure of the relevant indicators is completely described by a copula, including a nested Archimedean copula and its marginal tail behavior with peaks-over-threshold (POT) under the framework of extreme value theory.
• Second, in our pricing setting, the principal-based coupon paid out is closely related to the inter-annual and intra-annual variations of both the severity and frequency of the catastrophe risks with operational flexibility via the utilization of functions (see details in Section 2.1).
• Third, instead of modelling the occurrence process of disasters, we employ time series models to describe the heterogeneity of catastrophic events up to the maturity date of the CAT bonds. Note that the yearly paid coupon scenario facilitates the discretization of the counting process of catastrophe events, and time series including auto-regressive moving average models (ARMA) are helpful for modelling the time dependence of the quantities involved.
• Finally, we conduct an empirical analysis of pricing a rainfall CAT bond using comprehensive trigger indicators of the affected population, crop affected area, and direct economic losses in China during 2006–2020. Our pricing mechanism is well illustrated with both annual maxima and annual average retention levels (cf. Equations (10) and (11)). Monte Carlo simulations are carried out further to analyze the sensitivity of the rainstorm CAT bond price to trigger attachment levels, maturity date, catastrophe intensity, and numbers of trigger indicators (see details in Section 3).

The considerable CAT bond pricing mechanism is well-illustrated by pricing the flood CAT bond in China. The empirical analysis enables us to draw the following conclusions. First, the price of CAT bonds is negatively correlated with the maturity period, indicating that long-term bonds have lower prices due to the potential for more severe and frequent disasters, resulting in larger coupon and principal losses. Medium-term bonds tend to have relatively lower prices. Second, the time of bond purchase also affects the bond price, with earlier purchases leading to higher prices as they are associated with lower risks in previous years. Third, the trigger levels and intensity measures of the major rainstorms play significant roles in determining the bond prices. Higher trigger levels lead to increased coupon and principal retention, resulting in higher bond prices, while higher intensity measures of rainstorms lead to decreased bond prices due to reduced retained amounts. Finally, we found that CAT bonds triggered by multiple events are more attractive. Those CAT bonds triggered by more events have lower prices since they are more likely to have a smaller concurrent trigger probability with a larger principal retention.

The remainder of this paper is organized as follows. Section 2 presents the CAT bond pricing models with an innovative trigger mechanism. Section 3 gives an empirical analysis of flood CAT bond pricing in China and a numerical analysis of pricing sensitivity. We conclude this paper in Section 4 with extensional discussion. The relevant methodology and proofs are deferred to Appendix A.

2. Model Formulation

2.1. CAT Bond Pricing Model

This paper considers the pricing of catastrophe bonds due to a single disaster, where both the coupon and the principal are at risk in the case of a serious disaster. We refer to the main idea in Chao and Zou (2018) to consider a coupon-paying CAT bond triggered by m dependent catastrophe indicators $x_1,\dots ,x_m$. The investors may receive a portion of the coupon at the end of each year and a portion of the principal back at the maturity date. These proportions are determined by the accumulated excesses of the trigger indicators $x_i$ over their attachment levels $u_i$. Different from Chao and Zou (2018), we incorporate a triple of quantities $({\alpha }_t,{\beta }_t,{\gamma }_t)$ to determine the step-wise coupon and principal triggered. The triple plays an important role in the pricing of the CAT bond, which depends on the stressful indicators vector $\mathit{x}=(x_1,\dots ,x_m)$ and the counting process of the disaster $\left\{N\right(t),t\ge 0\}$, i.e., an integer-valued, non-negative, and non-decreasing stochastic process. We will develop a pricing mechanism reflecting both exceeding magnitude of the trigger indicators and its occurrence of the potential disaster below. We list all notations with explanations involved in

Table 1. Notation and its description involved in the pricing formula.

Notation	Description
F	Principal
$C_t$	Coupon paid in year $t,t=1,\dots ,T$
R	Coupon rate
$N\left(t\right)$	Number of disasters up to year t
$N_t=N\left(t\right)-N(t-1)$	Number of disasters in year t
$x_{ij}$	The jth observation of the ith trigger indicator up to year t, $j=1,\dots ,N\left(t\right)$
$u_i$	Attachment level of the ith trigger indicator, $i=1,\dots ,m$

.

Supposing that the observations of the m-dimensional indicator vector $(x_{1j},\dots ,x_{mj}),$$j=N(t-1)+1,\dots ,N\left(t\right)$ are independent of the counting process $\left\{N\right(t),t\ge 0\}$(Chao and Zou 2018; Goda 2021; Ibrahim et al. 2022), we define the overall catastrophe risk severity in year t as below:

$$\begin{array}{c}\hfill {\alpha }_t=f\left(s_{N(t-1)+1},\dots ,s_{N\left(t\right)}\right),s_j=\prod _{i=1}^m\left[1-\frac{{\left(x_{ij}-u_i\right)}_{+}}{x_{ij}}\right],\end{array}$$

(1)

where $x_{+}=max(x,0),f:{[0,1]}^{\mathbb{N}}\mapsto [0,1]$, a component-wise, non-decreasing, non-negative, and right-continuous function defined on a filtered physical probability space. We make f a general function so as to allow for different designs of CAT bonds. In general, we may suppose the coupon retention ${\alpha }_t$ to be the yearly average of the cumulative non-excess proportion caused by each disaster or its maxima as below.

(i)

The average coupon retention proportion due to the disaster in year t is modeled by

$${\alpha }_t=\frac{1}{N\left(t\right)-N(t-1)}\sum _{j=N(t-1)+1}^{N\left(t\right)}\prod _{i=1}^m\left[1-\frac{{\left(x_{ij}-u_i\right)}_{+}}{x_{ij}}\right].$$

(ii)

The maximum coupon retention proportion due to the disaster in year t is modeled by

$${\alpha }_t=\underset{N(t-1)+1\le j\le N\left(t\right)}{max}\prod _{i=1}^m\left[1-\frac{{\left(x_{ij}-u_i\right)}_{+}}{x_{ij}}\right].$$

Remark 1. 

Note that each $s_j\in (0,1]$. The case with $s_j=1$ shows that in the jth disaster, all trigger indicators are below the attachment levels. The smaller $s_j$ is, the more likely indicator $x_j=(x_{1j},\dots ,x_{mj})$ is to be far larger than its attachment level $u=(u_1,\dots ,u_m)$. It follows by the component-wise non-decreasing property of the function f that the ${\alpha }_t\in (0,1]$ is appropriate to quantify the coupon retention proportion.

To further trigger partial principal caused by multiple stressful indicators, i.e., more than one indicator being over its attachment level, we define the following two new indices ${\beta }_t$ and ${\gamma }_t$. Similar to $s_j$ and the function f in Equation (1), we define

$$\begin{array}{ccc}\hfill {\beta }_t& =& g\left(s_{N(t-1)+1}^{*},\dots ,s_{N\left(t\right)}^{*}\right),s_j^{*}=\prod _{1\le i_1<i_2\le m}\left[1-\frac{{\left(x_{i_{1}j}-u_{i_1}\right)}_{+}}{x_{i_{1}j}}·\frac{{\left(x_{i_{2}j}-u_{i_2}\right)}_{+}}{x_{i_{2}j}}\right],\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill {\gamma }_t& =& h\left(s_{N(t-1)+1}^{**},\dots ,s_{N\left(t\right)}^{**}\right),s_j^{**}=\prod _{1\le i_1<i_2<i_3\le m}\left[1-\frac{{\left(x_{i_{1}j}-u_{i_1}\right)}_{+}}{x_{i_{1}j}}·\frac{{\left(x_{i_{2}j}-u_{i_2}\right)}_{+}}{x_{i_{2}j}}·\frac{{\left(x_{i_{3}j}-u_{i_3}\right)}_{+}}{x_{i_{3}j}}\right],\hfill \end{array}$$

(3)

where g and h are two functions that are component-wise, non-decreasing, non-negative, and right-continuous functions from ${[0,1]}^{\mathbb{N}}$ to $[0,1]$. The following proposition gives the properties of $(s_j,s_j^{*},s_j^{**})$ and $({\alpha }_t,{\beta }_t,{\gamma }_t)$.

Proposition 1. 

Let $(s_j,s_j^{*},s_j^{**})$ be defined by Equations (1)–(3) and $({\alpha }_t,{\beta }_t,{\gamma }_t)$ be in $[0,1]$ obtained by the three component-wise, non-decreasing, non-negative, and right-continuous functions $(f,g,h)$ from ${[0,1]}^{\mathbb{N}}$ to $[0,1]$.

(i) 

It follows that $s_j,s_j^{*},s_j^{**}$ range over $[0,1]$. The case with $s_j^{*}<1$ (or $s_j^{**}<1$) implies that in the jth disaster, among all m trigger indicators, there are at least two (or three) indicators above its attachment levels simultaneously. Similar to $s_j$, the values of $s_j^{*}$ and $s_j^{**}$ quantify the concurrent exceeding magnitude over its thresholds.

(ii) 

Quantitatively, we have $0<s_j\le s_j^{*}\le s_j^{**}\le 1$. It follows further that

$$0<{\alpha }_t\le {\beta }_t\le {\gamma }_t\le 1,$$

when the three functions $f,g$, and h are taken as the same.

(iii) 

It follows by the component-wise non-decreasing monotonicity of the functions $f,g$, and h that the quantities ${\alpha }_t,{\beta }_t,{\gamma }_t$ are negatively associated with both the severity and the frequency of the disaster in the given year t.

Proposition 1 is key for our new pricing mechanism. The triple of resulting quantities $({\alpha }_t,{\beta }_t,{\gamma }_t)$ gives us the insight to trigger a partial coupon and principal when more indicators are triggered simultaneously. The proof of Proposition 1 is deferred to Appendix A. We will design below the cashflow of a CAT bond based on Proposition 1. Suppose that the investor buys a CAT bond with a face value of F in year t and maturity in year T. Denote by $C_{t,s}$ the coupon paid in year s and $F_{t,T}$ the redemption value in year T. We define

$$\left\{\begin{array}{c}{C}_{t,s}=\left\{\begin{array}{cc}{\alpha }_{t+1}·C_0,\hfill & s=t+1,\hfill \\ {\alpha }_s·\prod _{k=t+1}^{s-1}\left(\frac{{\beta }_k}{2}{\mathbb{I}}_{\left\{{\gamma }_k<1\right\}}+{\beta }_k{\mathbb{I}}_{\left\{{\gamma }_k\ge 1\right\}}\right)·C_0,\hfill & s=t+2,\dots ,T,\hfill \end{array}\right.\hfill \\ F_{t,T}=\prod _{k=t+1}^T\left(\frac{{\beta }_k}{2}{\mathbb{I}}_{\left\{{\gamma }_k<1\right\}}+{\beta }_k{\mathbb{I}}_{\left\{{\gamma }_k\ge 1\right\}}\right)·F,\hfill \end{array}\right.$$

(4)

where $({\alpha }_t,{\beta }_t,{\gamma }_t)$ are defined by Equations (1)–(3) and $C_0=F·R$, the coupon paid according to a fixed coupon rate R for a bond with the face value of F. We can see from Equation (4) that the coupon paid will be ${\alpha }_t$ times the initial coupon if there is at most one indicator triggered; otherwise, the principal will be wiped out and the coupon will be reduced accordingly. In this scenario, there are two cases: the first is that exactly two indicators are triggered simultaneously (i.e., ${\beta }_t<{\gamma }_t=1$), and the second is that at least three indicators are triggered simultaneously, i.e., ${\gamma }_t<1$. We differentiate these two cases with principal wipe-out by ${\beta }_t$ and ${\beta }_t/2$, respectively. This measure safeguards the initial investment capital of investors by ensuring that a portion of the principal remains preserved (Tang and Yuan 2019), thus increasing the investor’s interest in the investment of the CAT bond.

Remark 2. 

(a) 

We see that the coupon retained in year t is closely related to inter-annual variation of the catastrophe risks via the accumulative principal ${\prod }_{k<t}\left(\frac{{\beta }_k}{2}{\mathbb{I}}_{\left\{{\gamma }_k<1\right\}}+{\beta }_k{\mathbb{I}}_{\left\{{\gamma }_k\ge 1\right\}}\right)$, as well as the intra-annual risk via ${\alpha }_t$ with operational flexibility via the induction of functions $f,g$, and h.

(b) 

Our trigger mechanism in Equation (4) with the hierarchical proportional coupon and principal paid out may attract more investors in comparison with the hybrid trigger mechanism given by Wei et al. (2022), since therein the current and future coupons will be paid out once one of the indicators is triggered and the principal at maturity will be completely wiped out once both indicators are triggered simultaneously.

Finally, the price of the CAT bond, as the present value of all future cashflows, denoted by $P_t$, is given by

$$P_t=\sum _{s=t+1}^T\mathbb{E}\left\{C_{t,s}\right\}p(t,s)+\mathbb{E}\left\{F_{t,T}\right\}·p(t,T),$$

(5)

where $p(t,s)$ represents the discount factor at time t of the zero-coupon bond with a redemption value of 1 at time $s,t<s\le T$, given by

$$p(t,s)=\mathbb{E}\left\{exp\left(-{\int }_t^{s}r\left(u\right)\mathrm{d}u\right)\right\}$$

with interest rate process $\left\{r\right(t),t\ge 0\}$. In the following, we consider a CIR model for the interest rate risk, which will lead to an explicit expression of the discount factor.

2.2. Interest Rate Model

In our study, we adopt the Cox–Ingersoll–Ross (CIR) model to describe the continuous-time interest rate process $\left\{r\right(t),t\ge 0\}$ (Cox et al. 1985). The CIR model offers specific advantages over the Vasicek interest rate model, such as its simplicity, ease of handling, mean-reverting behavior, and avoidance of negative interest rates (Mistry and Lombardi 2022; Vasicek 1977). It satisfies the following stochastic differential equation under the risk-neutral probability measure:

$$dr\left(t\right)=k\left(\theta -r\left(t\right)\right)dt+\epsilon \sqrt{r\left(t\right)}d{W}_t,$$

(6)

where $k>0$ is the speed of mean reverting, $\epsilon >0$ is the volatility parameter, $\theta >0$ is the long-run mean of the interest rate, and {$W_t,t\ge 0\}$ denotes the standard Wiener process. Consequently, the discount factor $p(t,T)$ is given by

$$p(t,T)=A(t,T){\mathrm{e}}^{-B(t,T)r\left(t\right)},$$

where

$$\begin{array}{cc}\hfill A(t,T)& ={\left[\frac{2\eta {\mathrm{e}}^{(\kappa +\eta )(T-t)/2}}{(\kappa +\eta )\left({\mathrm{e}}^{\eta (T-t)}-1\right)+2\eta }\right]}^{2\kappa \theta /{\epsilon }^2},\hfill \\ \hfill B(t,T)& =\frac{2\left({\mathrm{e}}^{\eta (T-t)}-1\right)}{(\kappa +\eta )\left({\mathrm{e}}^{\eta (T-t)}-1\right)+2\eta },\hfill \\ \hfill \eta & =\sqrt{{\kappa }^2+2{\epsilon }^2}.\hfill \end{array}$$

2.3. POT Model

The extreme value theory plays an important role in analysing statistical patterns of extreme-value events. There are two typical approaches to extract the extreme samples: the block maximum (BM) method and the peaks-over-threshold (POT) model. As the POT model can make full use of the extreme-value data in comparison with the BM model, it is widely used in the fields of insurance, hydrology, and finance (Ma et al. 2021).

Suppose that $X_1,X_2,\dots ,X_n,\dots$ is a random sample from parent $X\sim F\left(x\right)$, i.e., the $X_i$ values are independent and identically distributed with a common distribution function (df) $F\left(x\right)$. Given a threshold u, the distribution $F_u\left(y\right)$ of the excess $Y_{\left[u\right]}=X-u|X>u$ is thus given by

$$F_u\left(y\right)=\mathbb{P}\left\{X-u\le y\mid X>u\right\}=\frac{\mathbb{P}\left\{u<X\le y+u\right\}}{\mathbb{P}\left\{X>u\right\}}=\frac{F(y+u)-F\left(u\right)}{1-F\left(u\right)},y>0.$$

For sufficiently high thresholds, Pickands (1975) pointed out that the distribution of the threshold-excess threshold $Y_{\left[u\right]}$ can be approximated by a generalized Pareto (GP) distribution $G_{\xi ,\sigma }(·)$, i.e., for the right endpoint $x^{*}=sup\{x\in \mathbb{R}:F\left(x\right)<1\}$

$$\begin{array}{c}{\displaystyle \underset{u\to x^{*}}{lim}\underset{0\le y\le x^{*}-u}{sup}}\left|F_u\left(y\right)-G_{\xi ,\sigma }\left(y\right)\right|=0.\end{array}$$

Therefore, the tail distribution function $\overline{F}\left(x\right)=1-F\left(x\right)$ of X can be approximated by

$$\overline{F}\left(x\right)=\overline{F}\left(u\right){\overline{F}}_u(x-u)\approx \overline{F}\left(u\right){\overline{G}}_{\xi ,\sigma }(x-u),x>u.$$

(7)

Here $\xi \in \mathbb{R}$ and $\sigma >0$ are the shape and scale parameters of the GP distribution $G_{\xi ,\sigma }\left(y\right)=1-{[1+\xi y/\sigma ]}^{-1/\xi },y>0$. In practice, the exceedance probability $\overline{F}\left(x\right)$ gives the insight into the potential risk. Its estimate can be obtained through the extrapolation approach via Equation (7): we obtain the approximate tail probability of the GP model using the maximum likelihood estimation of $\xi ,\sigma$ based on the excesses $(x_{\left(i\right)}-u)$ with $x_{\left(1\right)}\ge \cdots \ge x_{\left(n_u\right)}$ exceeding the threshold u and the estimate of $\overline{F}\left(u\right)$ as $n_u/n$. Theoretically, the threshold u can be determined by minimizing the mean square error of the Hill estimator of $\xi$, balancing the model bias and variance. A common graphical approach to the determination of the threshold is to check both the linearity of the empirical mean excess function

$$\begin{array}{c}\hfill e_n\left(u\right)=\frac{1}{n_u}\sum _{i=1}^{n_u}\left(x_{\left(i\right)}-u\right)\end{array}$$

(8)

and its derived stable estimates of both scale and shape parameters.

2.4. Copula

The concept of a copula, proposed by Abe (1973), serves as a tool for describing the dependence structure among the marginal variables. Namely, for a joint distribution $G(x_1,\dots ,x_m)$ with marginal df $G_i$ for the ith component, the copula C is thus determined by

$$C(u_1,\dots ,u_m)=G(G_1^{-1}\left(u_1\right),\dots ,G_m^{-1}\left(u_m\right)),(u_1,\dots ,u_m)\in {[0,1]}^m.$$

Table 2. Common Archimedean copulas.

Copula	Generator	$\mathit{C}({\mathit{u}}_1,\mathit{\dots },{\mathit{u}}_{\mathit{m}};\mathit{\theta })$	Parameters
Gumbel	${(-lnt)}^{\theta }$	$exp\left(-{\left({\sum }_{i=1}^m{\left(-ln{u}_i\right)}^{\theta }\right)}^{1/\theta }\right)$	$\theta \ge 1$
Clayton	$\frac{t^{-\theta }-1}{\theta }$	${\left({\sum }_{i=1}^m{u}_i^{-\theta }-m+1\right)}^{-1/\theta }$	$\theta \ge 0$
Frank	$-ln\frac{{\mathrm{e}}^{-\theta t}-1}{{\mathrm{e}}^{-\theta }-1}$	$-\frac{1}{\theta }ln\left(1+{\prod }_{i=1}^m\frac{{\mathrm{e}}^{-\theta u_i}-1}{{\left({\mathrm{e}}^{-\theta }-1\right)}^{m-1}}\right)$	$\theta \ne 0$

lists common Archimedean copulas, including the Clayton, Gumbel, and Frank copulas, with a convex and decreasing generator $\varphi \left(t\right):[0,1]\mapsto (0,\infty )$ satisfying $\varphi \left(1\right)=0$. Its tail dependence is controlled by the parameter $\theta$ and its copula is given by

$$C\left(u_1,\dots ,u_m\right)={\varphi }^{-1}\left(\sum _{i=1}^m\varphi \left(u_i\right)\right).$$

Given the analytic tractability of Archimedean copulas, they are widely applied in insurance, finance, hydrology, survival analysis, etc. In this paper, we consider hierarchical (or nested) Archimedean copulas representing the different dependence among the components (Hofert 2010). Namely, the nested Archimedean copula is of the form

$$\begin{array}{c}\hfill C(u_1,\dots ,u_m)=C_{outer}\left(C_{inner}\left(u_1,\dots ,u_k;{\theta }_1\right),u_{k+1},\dots ,u_m;{\theta }_2\right),\end{array}$$

(9)

where the inner and outer copulas could be one of the three Archimedean copulas in Table 2. The tail dependence could be measured by the parameters ${\theta }_1$ and ${\theta }_2$. As shown in Charpentier and Segers (2009), the larger $\theta$ involved in the Gumbel, Clayton, or Frank copulas indicates a stronger dependence among these variables.

Combining the above pricing model and statistical model, this paper derives a rainstorm catastrophe bond. An empirical analysis together with a numerical analysis is given in Section 3.

3. Empirical Analysis and Numerical Analysis

To illustrate our pricing mechanism, this section develops the pricing of a severe rainstorm in China based on all recorded 245 major rainstorms in China during 2006–2020 and three main hazard indicators from the China Statistical Yearbook of Natural Disasters. These indicators cover all the disaster indicators, namely, the affected population (AP) accumulating the number of deaths, missing, emergency re-locations, and people with drinking water damaged; the crop affected area (CAA), reflecting both the cropland flooding area and its harvest-affected area; and the direct economic loss (DEL) adjusted with the Consumer Price Index (CPI) in 2020, summing the damaged losses of collapsed/damaged houses and other properties multiplied by its damaged ratio.

In what follows, we will first present the descriptive statistics of our triple of trigger indicators $(X_1,X_2,X_3)$ representing AP, CAA, and DEL in Section 3.1. The non-normality of these indicators motivates us to split the full range of the data into bulk and tail parts using classical parametric models and a peaks-over-threshold (POT) approach in Section 3.2, followed by the investigation of the joint distribution functions for disaster indicator data using a nested Archimedean copula in Section 3.3. The ARMA model is used to predict the disaster intensity over the next 3 years in Section 3.4. Monte Carlo simulation of the CAT bond price is carried out in Section 3.5. Here, we take the retention functions of $f,h$, and g in Equations (1)∼(3) as the same as the average or maximum function unless stated otherwise. That is,

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\alpha }_t\hfill & =\frac{1}{N\left(t\right)-N(t-1)}\sum _{j=N(t-1)+1}^{N\left(t\right)}\prod _{i=1}^m\left[1-\frac{{\left(x_{ij}-u_i\right)}_{+}}{x_{ij}}\right].\hfill \\ {\beta }_t\hfill & =\frac{1}{N\left(t\right)-N(t-1)}\sum _{j=N(t-1)+1}^{N\left(t\right)}\prod _{1\le i_1<i_2\le m}\left[1-\frac{{\left(x_{i_{1}j}-u_{i_1}\right)}_{+}}{x_{i_{1}j}}·\frac{{\left(x_{i_{2}j}-u_{i_2}\right)}_{+}}{x_{i_{2}j}}\right],\hfill \\ {\gamma }_t\hfill & =\frac{1}{N\left(t\right)-N(t-1)}\sum _{j=N(t-1)+1}^{N\left(t\right)}\prod _{1\le i_1<i_2<i_3\le m}\left[1-\frac{{\left(x_{i_{1}j}-u_{i_1}\right)}_{+}}{x_{i_{1}j}}·\frac{{\left(x_{i_{2}j}-u_{i_2}\right)}_{+}}{x_{i_{2}j}}·\frac{{\left(x_{i_{3}j}-u_{i_3}\right)}_{+}}{x_{i_{3}j}}\right]\hfill \end{array}\right.\end{array}$$

(10)

and

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\alpha }_t\hfill & =\underset{j=N(t-1)+1}{\overset{N\left(t\right)}{max}}\prod _{i=1}^m\left[1-\frac{{\left(x_{ij}-u_i\right)}_{+}}{x_{ij}}\right].\hfill \\ {\beta }_t\hfill & =\underset{j=N(t-1)+1}{\overset{N\left(t\right)}{max}}\prod _{1\le i_1<i_2\le m}\left[1-\frac{{\left(x_{i_{1}j}-u_{i_1}\right)}_{+}}{x_{i_{1}j}}·\frac{{\left(x_{i_{2}j}-u_{i_2}\right)}_{+}}{x_{i_{2}j}}\right],\hfill \\ {\gamma }_t\hfill & =\underset{j=N(t-1)+1}{\overset{N\left(t\right)}{max}}\prod _{1\le i_1<i_2<i_3\le m}\left[1-\frac{{\left(x_{i_{1}j}-u_{i_1}\right)}_{+}}{x_{i_{1}j}}·\frac{{\left(x_{i_{2}j}-u_{i_2}\right)}_{+}}{x_{i_{2}j}}·\frac{{\left(x_{i_{3}j}-u_{i_3}\right)}_{+}}{x_{i_{3}j}}\right].\hfill \end{array}\right.\end{array}$$

(11)

3.1. Descriptive Analysis of Trigger Indicators

We can see from

Table 3. Descriptive statistics of affected population (AP) in millions, crop affected area (CAA) in million hectares, and direct economic loss (DEL) in billion yuan.

Trigger Indicators	Maximum	Minimum	Mean	Median	Skewness	Kurtosis
AP	1510.000	1.600	187.880	116.750	2.771	12.678
CAA	135.200	0.005	12.492	7.000	3.689	20.458
DEL	420.816	1.011	21.411	8.083	5.154	37.823

that the averages of all three main indicators are far larger than their respective medians. Moreover, both the skewness and the kurtosis of these indicators, which are far larger than 0 and 3, respectively, indicate that the right-skewed pattern of these indicators have possible heavy tails. This is graphically consistent with the exponential QQ plots in Figure 2 with a downward convex deviation from the theoretical straight line.

3.2. POT-Based Tail Analysis of Trigger Indicators

In order to determine the threshold levels $(u_1,u_2,u_3)$ of the triple of indicators $(X_1,X_2,X_3)$, we examine the mean residual life plots in Figure 3 (see Equation (8) for details). We determine first the threshold $u_1$ for the affected population $X_1$. The possible range of the threshold is detected to be in the interval of $(150,200)$ according to the sample mean residual plot with a linear trend (see Figure 3a), and then we investigate the stability of the estimations of the scale and shape parameters involved in the generalized Pareto model of the threshold excesses of AP in Figure 4a,b. Consequently, the attachment level $u_1=160$ is determined for the affected population (AP). Similar arguments are given to the attachment levels of $u_2=12$ for the crop affected area (CAA) and $u_3=15$ for the direct economic loss (DEL).

We model the probability distribution of our trigger indicators using the full range of the data, with sufficient flexibility for separate control over the bulk and tail features. Different from Chao and Zou (2018), we consider the Beta-GP models for the scaled non-exceedances and threshold excesses, namely,

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{X}_i-u_i|X_i>u_i\sim G_{}\left(y\right)=1-{\left(1+\frac{\xi y}{\sigma }\right)}^{-1/\xi },y>0,\hfill & \mathrm{for}\mathrm{threshold}\mathrm{excesses},\hfill \\ X_i^{*}=\frac{X_i-m_i}{u_i-m_i}\sim Bet{a}_{a,b}\left(y\right)=\frac{1}{B(a,b)}{y}^{a-1}{(1-y)}^{b-1},0<y<1,\hfill & \mathrm{for}\mathrm{non-exceedances},\hfill \end{array}\right.\end{array}$$

(12)

where $(m_1,m_2,m_3)=(1.6,0.005,1.011)$ are the sample minima of the trigger indicators of AP, CAA, and DEL, respectively, given in

Table 4. Estimates of the scale $\left(\sigma \right)$ and shape $(\xi$) parameters in the GP models of the threshold excesses and the two parameters $\alpha$ and $\beta$ in the Beta model $Bet{a}_{\alpha ,\beta }$ for the scaled bulk data below the threshold. The p-value based on a Chi-square test confirms that the Beta-GP model fits the data well.

	GP					Beta
	Scale ($\mathit{\sigma }$)	95% CI	Shape ($\mathit{\xi }$)	95% CI	p-Value	a	b	p-Value
AP	173.369	(118.21, 226.38)	0.197	(−0.04, 0.44)	0.247	1.016	1.345	0.246
CAA	11.771	(7.41, 16.13)	0.341	(0.03, 0.64)	0.283	1.076	1.687	0.299
DEL	23.538	(13.96, 33.19)	0.492	(0.13, 0.84)	0.250	0.582	1.137	0.253

, and $(u_1,u_2,u_3)=(160,12,15)$ are given by the mean residual plots and parameter stability plots. All maximum likelihood estimations of the parameters involved in the Beta-GP models are shown in Table 4. We see that both CAA and DEL possess heavy tails with 95% confidence intervals of the shape parameters of $(0.03,0.64)$ and $(0.13,0.84)$, respectively, while the estimated shape parameter for AP excesses is 0.197, showing a certain power-decaying tail. Moreover, we examine the model goodness-of-fit using a Chi-square test and obtain that all the p-values are larger than 0.246, confirming that the Beta-GP model agrees with the observed trigger indicators. Intuitively, Figure 5 illustrates the appropriateness of the GP model of the threshold exceedances for each trigger indicator, since both the PP plots and QQ plots show that almost all points are around the straight line.

To conclude, the fitted distributions $F_i,i=1,2,3$ corresponding to the trigger indicators of the affected population (AP), the crop affected area (CAA), and the direct economic loss (DEL) are given below

$$\begin{array}{c}\hfill F_i\left(x\right)=\left\{\begin{array}{cc}1-\frac{n_{u_i}}{n}{\overline{G}}_{{\xi }_i,{\sigma }_i}(x-u_i),\hfill & x>u_i,\hfill \\ \frac{1}{B(a_i,b_i)}{\int }_0^{(x-m_i)/(u_i-m_i)}{t}^{a_i-1}{(1-t)}^{b_i-1}dt,\hfill & x\le u_i,\hfill \end{array}\right.\end{array}$$

(13)

where $n=245$ and the threshold $(u_1,u_2,u_3)$, the excess numbers $(n_{u_1},n_{u_2},n_{u_3})$, the sample minina $(m_1,m_2,m_3)$, and other parameters for the Beta-GP model are given by Equation (12) and Table 4.

Given the distributions of each trigger indicator in Equation (13), it remains to discuss the dependence via nested Archimedean copula in the following section in order to illustrate our multiple-event-triggered pricing mechanism in Section 3.5.

3.3. Dependence Anaylsis of Trigger Indicators Based on Nested Archimedean Copula

We shall investigate first the non-exchangeable dependence among the trigger indicators through the Spearman rank correlation $\rho$. The Spearman $\rho$ for (AP, CAA), (AP, DEL), and (CAA, DEL) are 0.771, 0.554, and 0.515, respectively. We see that all three pairs demonstrate a certain degree of dependence, with a stronger dependence within (AP, CAA) than between them. This motivates us to describe the dependence using a nested Archimedean copula, with an inner copula for (AP, CAA) with Archimedean parameter ${\theta }_1$ and DEL placed into the outer Archimedean copula with parameter ${\theta }_2$ (see Equation (9)).

Next, to model the dependence structure using a nested Archimedean copula, we construct a uniformly distributed transformation of the raw data according to the marginal analysis in Section 3.2, i.e., a straightforward application of Equation (13) given by

$${\tilde{u}}_{ij}=F_i\left(x_{ij}\right),i=1,2,3,j=1,2,\dots ,245.$$

Finally, we examine the dependence structure using a nested Archimedean copula with both the inner and outer copula being one of the Gumbel, Clayton, and Frank copulas listed in Table 2. We suppose that

$$\begin{array}{c}\hfill ({\tilde{u}}_{1j},{\tilde{u}}_{2j},{\tilde{u}}_{3j})\stackrel{i.i.d.}{\sim }{C}_{outer}(C_{inner}({\tilde{u}}_{1j},{\tilde{u}}_{2j};{\theta }_1),{\tilde{u}}_{3j};{\theta }_2),j=1,2,\dots ,245.\end{array}$$

(14)

We see from

Table 5. Maximum likelihood estimates of the parameters involved in the inner and outer copulas in Equation (9). The Kolmogorov–Smirnov (KS) test and its p-value indicates that the best model is the nested Gumbel copula.

Inner and Outer Copula	${\mathit{\theta }}_1$ (Inner)	${\mathit{\theta }}_2$ (Outer)	KS Test	p-Value
Gumbel	44.68	22.80	0.06	0.12
Clayton	9.58	4.35	0.12	0.35
Frank	176.34	87.60	0.04	0.61

that the maximum likelihood estimates of the inner parameter ${\theta }_1$ are all larger than that of the outer parameter ${\theta }_2$, agreeing with the stronger dependence within the inner variables than between them. Moreover, we will select the most suitable copula based on the Kolmogorov–Smirnov (KS) test. We see that all nested Archimedean copulas fit the data well and the best nested Gumbel copula is determined with the minimal KS test value of 0.04 and the maximum p-value of 0.61.

3.4. Modelling of Annual Frequency of Rainstorms in China

Note that our pricing mechanism is of discrete form. It follows from Equation (4) that the proportion of the coupon and principal paid out depends not only on the severity of the disasters (thus, the trigger indicators) but also on the independent annual frequency of disasters. This was confirmed by our Pearson correlation analysis between the annual number of rainstorms and its resulting annual average of losses in terms of AP, CAA, and DEL, with values equal to $0.08,-0.09$, and $-0.1$, respectively. Due to the changing climate, the occurrence of major rainstorms becomes more and more frequent.The purpose of this section is to model the intensity of annual major rainstorms during 1986–2020 in China using an auto-regressive moving average (ARMA) model so as to forecast the frequency in the next three years, i.e., 2021–2023. The relevant data are from the China Statistical Yearbook of Natural Disasters.

The stationarity of the time series ${\mathrm{\Lambda }}_t$ is validated by the augmented Dickey–Fuller (ADF) test at a significance level of 0.05. Next, we determine the autoregressive order and moving-average order as $p=1,q=3$, respectively, by the position of the lag cutoff in the partial autocorrelation function (PACF) and ACF plots in Figure 6. Thus, the ARMA(1, 3) model is selected for modelling the annual intensity of major rainstorms in China and the maximum likelihood estimates of the parameters involved are obtained as $({\varphi }_1,{\theta }_1,{\theta }_2,{\theta }_3)=(0.816,0.268,0.240,-0.748)$.

The Jarque–Bera (JB) test with a p-value of 0.08 confirms the normality assumption of the residuals $e_k$, given in Equation (A1) in the ARMA model, and its i.i.d. assumption is verified by the Ljung–Box (LB) test with $p=0.502$. Before we apply this model for prediction, we carried out an F-test, which confirmed the goodness of fit with $p=0.003$. Consequently, the straightforward application of ARMA(1,3) to the intensity measures for 2021, 2022, and 2023 predicts $(41.86,41.56,39.39)$ a general upward trend in intensity for the next three years, as shown in Figure 7, which will be used in Section 3.5 below for the numerical analysis of pricing sensitivity.

3.5. Pricing of CAT Bond and Numerical Analysis

This section will focus on the CAT bond pricing based on our pricing mechanism in Equation (5) and the basic analyses of major rainstorms we conducted in Section 3.1, Section 3.2, Section 3.3 and Section 3.4. In the current pricing mechanism, we suppose mutual independence among interest rate, disaster indicators, and counting process (Chao and Zou 2018; Goda 2021; Ibrahim et al. 2022; Reshetar 2008). This facilitates Monte Carlo simulations in both the numerical analysis and the empirical analysis. In practice, there might be weak and ignorable dependence between the financial market and natural disaster risks, which attracts investors for the benefit of diversification of their investments. On the other hand, the dependence assumption between the frequency and severity of disasters is more likely to bring about model mis-specification. In what follows, we consider a T-year period CAT bond with a principal of $F=CNY100$ and a fixed coupon rate $R=3.5\%$. Note that the explicit expectation in Equation (5) is not available. We thus simulate the price of the CAT bond and estimate it as the average of the simulated price. The steps of the simulation are outlined as follows.

(1)

Generate T random numbers $N_1,N_2,\dots ,N_T$ independently from a Poisson distribution with intensity measure ${\mathrm{\Lambda }}_1,{\mathrm{\Lambda }}_2,\dots ,{\mathrm{\Lambda }}_T$, standing for the number of major rainstorms in year $1,2,\dots ,T$. Here, $({\mathrm{\Lambda }}_1,{\mathrm{\Lambda }}_2,\dots ,{\mathrm{\Lambda }}_T)$ can be provided by the ARMA(1, 3) model in Section 3.4.

(2)

Based on the marginal distribution in Equation (13) and the nested Frank copula in Equation (14) with parameters given in Table 4 and Table 5, we generate a sample of $(x_{1k},x_{2k},x_{3k})$ of size $N_j$ for year $j,1\le j\le T$ as

$$x_{ik}=F_i^{-1}\left({\tilde{u}}_{ik}\right),i=1,2,3,k=\sum _{l=1}^{j-1}{N}_l+1,\dots ,\sum _{l=1}^j{N}_l.$$

(3)

We suppose that the stochastic interest rate process $\left\{r\right(t),t\ge 0\}$ follows the CIR model stated in Equation (6), and the parameter estimates involved are obtained through maximum likelihood estimation based on one-year seven-day interest rate data for 2021, available at the Shanghai Interbank Offered Rate. We select the interest rate on 6 January 2021 as the initial interest rate, i.e.,

$$\begin{array}{c}\hfill \left\{\begin{array}{c}dr\left(t\right)=0.2\left(0.05-r\left(t\right)\right)dt+0.05\sqrt{r\left(t\right)}d{W}_t,\hfill \\ r\left(0\right)=2.962\%.\hfill \end{array}\right.\end{array}$$

(15)

We simulate the interest risk to give the discount factor $p(t,T)$.

(4)

Given a trigger level $(u_1,u_2,u_3)$, we calculate the price of the CAT bond based on the cashflows in Equation (4) and the interest risk model in Equation (15).

We repeat $m={10}^4$ all the simulation steps above and obtain the sample mean of the price $P_t$ of a T-year CAT bond bought in year $t=0,1,\dots ,T-1$. In the following, we will discuss the pricing sensitivity in the maturity period, the trigger level, and the trigger indicators.

Price sensitivity in maturity period.

Table 6. CAT bond pricing of major rainstorms in China purchased in year t and maturating in year $T=1,2,3$ based on Equation (5). Here the trigger level $(u_1,u_2,u_3)$ is taken as the sample 90% quantile of (AP, CAA, DEL). The bond prices $P_{t1}$ and $P_{t2}$ are given based on the cashflows with an average retention level in Equation (10) and max retention level in Equation (11).

	$\mathit{T}=1$	$\mathit{T}=2$		$\mathit{T}=3$
	$\mathit{t}=\mathbf{0}$	$\mathit{t}=\mathbf{0}$	$\mathit{t}=\mathbf{1}$	$\mathit{t}=\mathbf{0}$	$\mathit{t}=\mathbf{1}$	$\mathit{t}=\mathbf{2}$
Bond price $P_{t1}$	80.392	70.881	66.255	56.006	54.865	52.321
Bond price $P_{t2}$	83.259	73.211	68.450	58.812	57.649	53.984

shows a T-year CAT bond triggered by the triple of trigger indicators (AP, CAA, DEL) with trigger level $(u_1,u_2,u_3)$ being the sample 90% quantile. The bond price decreases at the maturity period T and further, the downward trend in the bond price gradually increases. As the bond is issued for a long maturity period, future disasters might be more severe and frequent, causing a larger magnitude of triggered events and thus a smaller proportion of the coupon and principal retained. Consequently, a lower price is obtained for a medium-term CAT bond. Meanwhile, for a given maturity period $T=2$ or 3, the CAT bond might be sold in different years during the bond issue period. We see that the price of a bond purchased in an earlier year is higher, as the potential risk of major rainstorms might be lower in the previous few years and the coupon and principal are more likely to remain. Therefore, the pricing mechanism is fairly attractive in terms of raising capital for the purpose of reinsurance. We see further that the CAT bond prices vary in coupon–principal retention functions. Clearly, the coupon and principal remaining based on the annual average retention level in Equation (10) are not larger than that based on the annual maximum retention level in Equation (11), which results in the lower price of the former. Hence, our functional mechanism of the retention levels provides sufficient flexibility to the CAT bond issuers.

Price sensitivity to trigger level and intensity measure of major rainstorms. An appropriate trigger level is important for balancing the benefits between the investors and the bond issuers by means of the CAT bond price. In general, both the trigger levels and the intensity measures are determined by the severity and frequency of the potential disasters (here rainstorms), which can cause a change in the CAT bond price via the wiped-off coupon and principal in Equation (4). Clearly, we see from Figure 8 that the bond prices are positively correlated with the trigger levels but negatively associated with the intensity measures. Indeed, as the trigger level increases, a bond is less likely to be triggered and thus the coupon and principal retention level will increase, leading to an increase in the CAT bond price. In contrast, if major rainstorms occur more frequently with a larger intensity measure, the triggered indicators may accumulate a small amount of coupon and principal retained, thus causing a decrease in bond prices.

Price sensitivity in the selection of trigger indicators. As mentioned before, the multiple-event-triggered CAT bond receives increasing attraction from both investors and CAT issuers. In

Table 7. Comparison of one-year CAT bond price with different selections of multiple-event triggers under our pricing mechanism in Equation (5). The bond prices $P_{t1}$ and $P_{t2}$ are given according to the annual average and annual maximum retention functions in Equation (10) and Equation (11), respectively.

Trigger Indicator	AP-CAA-DEL	AP-CAA	AP-DEL	CAA-DEL
Bond price $P_{t1}$	80.3922	48.7558	48.7561	48.7564
Bond price $P_{t2}$	83.2592	48.7869	48.7871	48.7873

, we compare the one-year bond prices with different pairs of trigger indicators under our pricing mechanism in Equation (5). Apparently, the three-event-triggered CAT bond is issued with a lower price than that of two-event-triggered ones. This is because the smaller trigger probability of concurrent trigger events leads to a lower expectation of discounted cashflows. Additionally, the principal might be half-returned in case all three indicators are triggered (recalling the utilization of ${\gamma }_t$).

4. Conclusions and Extensional Discussions

In this paper, a multiple-event-triggered CAT bond pricing model is designed and a pricing formula is derived based on the copula–POT model. Our pricing model is more flexible and of more practical significance in terms of possessing a dynamic association between the coupon and principal paid off, whose magnitude varies in both the severity and frequency of the potential disasters. This model may provide a certain reference for the pricing research and subsequent practical application of catastrophe bonds. Although the multi-event trigger mechanism involves more elements and more complex processes and requires higher technical capabilities, the multiple-event-triggered bonds receive increasing attraction due to their low moral hazard and trigger risk (Chao and Zou 2018; Ibrahim et al. 2022; Wei et al. 2022). Meanwhile, it turns out that the CAT bond price decreases in the bond maturity period in the simulation, and it is negatively correlated with the trigger level and positively correlated with the catastrophe intensity.

The pricing mechanism depends on the joint distribution of the target indicators caused by the disaster (Goda 2021). In real life, catastrophic risks involve multiple disasters over many regions. Thus, a forthcoming research direction is to design a regional CAT bond pricing model that incorporates spatio-temporal extremes and the interaction of simultaneous disasters. Furthermore, we may bring the insights from Eling and Wirfs (2019) and incorporate relevant covariates such as the economic development level and the vulnerability to real rainstorms for estimating the local parameter of the POT model, which thus enables a region-specific pricing mechanism in practice. In addition, to relax the independence assumption among the interest rate risk, disaster indicators, and counting process might be more realistic in the CAT bond pricing models, with potential complexity regarding the pricing simulations and statistical tests of the possible dependence structure. Furthermore, the specific assumptions should be carefully verified since otherwise it is more likely to result in a model mis-specification issue. In practice, there might be weak dependence between the financial market risk and disaster risk, which hopefully enables the investors to be protected from the diversification investment.