Earthquake Bond Pricing Model Involving the Inconstant Event Intensity and Maximum Strength

Abstract

Traditional insurance’s earthquake contingency costs are insufficient for earthquake funding due to extreme differences from actual losses. The earthquake bond (EB) links insurance to capital market bonds, enabling higher and more sustainable earthquake funding, but challenges persist in pricing EBs. This paper presents zero-coupon and coupon-paying EB pricing models involving the inconstant event intensity and maximum strength of extreme earthquakes under the risk-neutral pricing measure. Focusing on extreme earthquakes simplifies the modeling and data processing time compared to considering infinite earthquake frequency occurring over a continuous time interval. The intensity is accommodated using the inhomogeneous Poisson process, while the maximum strength is modeled using extreme value theory (EVT). Furthermore, we conducted model experiments and variable sensitivity analyses on EB prices using earthquake data from Indonesia’s National Disaster Management Authority from 2008 to 2021. The sensitivity analysis results show that choosing inconstant intensity rather than a constant one implies significant EB price differences, and the maximum strength distribution based on EVT matches the data distribution. The presented model and its experiments can guide EB issuers in setting EB prices. Then, the variable sensitivities to EB prices can be used by investors to choose EB according to their risk tolerance.

1. Introduction

Extreme earthquakes are rare events and affect the economy and human life [1,2,3]. An analysis of data from the International Disaster Database (accessed on 28 October 2022 from https://www.emdat.be) shows that, on average, a single extreme earthquake results in economic damage of USD 1,266,348,909, with the highest magnitude compared to other extreme disasters, as illustrated in Figure 1. This high loss burdens the budgets of earthquake-prone countries (countries whose territories lie atop tectonic plate lines) [4,5,6,7], so several of them used traditional insurance mechanisms to obtain contingency funds for mitigation from 1992 to 2000. However, this mechanism generally has shortcomings, whereby the provided contingency costs pale in comparison to the actual losses incurred [8,9,10]. For example, earthquake insurance in Turkey in 1992 and Japan in 1993 provided contingency costs of USD 10.8 million (1.44% of actual losses) [11] and USD 16 million (1.6% of actual losses) [12], respectively. Therefore, to overcome this shortcoming, in the last three decades, earthquake-prone countries have attempted to develop traditional earthquake insurance mechanisms to be more effective in providing post-earthquake contingency costs [13,14,15]. So far, this development has established an earthquake insurance mechanism connected to financial instruments in the capital market. Simply put, investors receive a portion or the entirety of the earthquake risk associated with a country through financial instruments traded on the capital market. This is known as earthquake-insurance-linked securities (EILS).

In EILS, bonds are the most successful financial instruments because these can raise significant funds quickly with moderate risk [16,17]. This is what is called an earthquake bond (EB). Several earthquake-prone countries have issued EB as a contingency cost source. The history of countries using EBs is given in Table 1. Table 1 shows that most of the EBs were issued by countries on the American continent. Only one was not—in the Philippines. The first EB was issued by Mexico in 2006 [18,19,20]. Then, Mexico’s steps were followed by Latin American countries, namely, Chile, Colombia, and Peru, in 2018 [21]. Finally, the Philippines was the last to publish it in 2019 [22]. The term of an EB is generally two or three years with a contingency fund of USD 75 to 500 million. This contingency fund amount is more significant than that found in traditional insurance mechanisms [23], such as those in Turkey and Japan, which were mentioned previously. This then reduces the possibility of disproportion between the two in the insured country.

Based on Table 1, we see that a few earthquake-prone countries still issue EBs. One of the reasons is that the pricing mechanism is still being studied, e.g., by Indonesia. Determining the fair price of EBs is a fundamental stage of EB issuance, which is complicated because it integrates financial and earthquake risk variables, such as loss, strength, or the intensity of the earthquake [24,25].

Several studies have been conducted regarding modeling catastrophe (CAT) bond prices in general. Jarrow [26] modeled CAT bond prices in closed form using a robust model to obtain the solution analytically. Then, Nowak and Romaniuk [27] designed a CAT bond pricing model assuming independence between CAT and interest rate risks using the Hull–White and Cox–Ingersoll–Ross (CIR) models. Then, Ma and Ma [28] designed a CAT bond pricing model with no closed-form solution using an inhomogeneous compound Poisson process. The solution of the model is sought using an approximation method proposed by Chaubey et al. [29]. Then, Liu et al. [30] designed a CAT bond pricing model considering credit risk factors using Jarrow and Turnbull’s model [31]. Ma et al. [32] modeled CAT bond pricing, whose catastrophe characteristics are accommodated by integrating the Black Derman Toy [33] model and a double compound Poisson process. Then, Nowak and Romaniuk [34] introduced a CAT bond pricing model using a diffusion process with jumps and a multi-factor CIR model. In their research, Chao and Zou [35] developed a CAT bond pricing model based on multiple trigger events with a loss and fatality index.

Instead of modeling CAT bonds for general types of disasters, other studies have explored earthquake bond pricing models, with several models being developed and applied to various regions. Romaniuk [14] introduced a simple EB price model and a stochastic iteration method as an alternative way to obtain the model solution. Zimbidis et al. [8] introduced an EB price model with a discrete period using the extreme value theory (EVT) approach, the Vasicek stochastic interest rate model, and geometric Brownian motion. Then, the model is applied to earthquake data in Greece. Shao et al. [36] added earthquake depth factors and real interest rates to the Zimbidis et al. [8] model, which was applied to earthquake data in the United States. Gunardi and Setiawan [37] used the block maxima method and the CIR model to categorize EB prices into three types: zero coupon, coupon at risk, and principal and coupon at risk. The model is applied to earthquake data in Indonesia. Tang and Yuan [13] combined distorted and risk-neutral probability measures to design earthquake and interest rate risks. This model was then used for experiments estimating EB prices in Mexico and the United States. Then, Hofer et al. [38] introduced a municipal EB pricing model involving spatial weights between provinces in Italy, while Wei et al. [10] designed an EB price model based on a double trigger index (earthquake loss and strength) using the EVT and copula approaches. Mistry and Lombardi [39] developed the approach of Hofer et al. [38] with a constant coupon and a constant earthquake intensity rate, while Kang et al. [40] designed an EB pricing model considering the level of liabilities described through the asset liability management model. Aghdam et al. [41] introduced an EB pricing model using a Chebyshev-based spectral method. Anggraeni et al. [42] designed a municipal EB pricing model using EVT and the K-Means clustering method.

In reality, extreme earthquake intensity is inconstant across different time intervals. However, in the eleven related studies on EB price modeling explained in the previous paragraph, the risk of earthquake strength is still formulated based on a constant earthquake intensity (in fact, Hofer et al. [38] used one that is inconstant; however, it was used to describe the risk of loss, not strength). Therefore, previous models have not been able to accommodate this inconstant intensity, and this represents a gap that substantiates the novelty of this study.

Based on the research gaps presented, this research aims to design an EB pricing model that involves the inconstant event intensity and maximum strength of extreme earthquakes under the risk-neutral pricing measure. The inconstant intensity is measured based on only the finite frequency of extreme earthquake events. It is done to simplify the modeling stage and data processing time compared to considering all earthquake events that continuously occur throughout time. The inconstant intensity is formulated by an inhomogeneous Poisson process, while the maximum strength is modeled using extreme value theory (EVT). Then, the model is used in experiments using earthquake data in Indonesia, the country with the second highest frequency of earthquakes in the world. Finally, the variable sensitivities to the EB price involved in the model are also analyzed. The maximum strength risk model for extreme earthquakes can be used by practitioners, disaster management agencies, and geological agencies to measure the risk of earthquake severity based on its strength in a region. Then, the EB price models can be used by the issuer in setting reasonable EB prices. Then, the experiments conducted can guide us in how to use the models designed in this research. Investors can use the sensitivity of the variables to EB prices to choose an EB that suits their risk tolerance.

2. A Brief Explanation of Earthquake Bond

An earthquake bond is a type of earthquake insurance-linked security (ILS), where a sponsor, such as an insurer, reinsurer, or government, partially or fully offloads earthquake risk onto investors in exchange for elevated yields. This risk transfer is done by establishing a special purpose vehicle (SPV), which protects the sponsor and raises funds by issuing bonds to investors with a term of one to five years [43]. The SPV acquires bond and premium proceeds from investors and sponsors, respectively. Then, the SPV allocates them to secure financial instruments, such as Treasury bonds, and deposits the proceeds into a trust account. The proceeds are typically exchanged for variable returns linked to the London Interbank Offered Rate (LIBOR) provided by a reputable counterparty engaged in swap transactions. The purpose of the swap is to protect the sponsor and investors against potential risks associated with interest rate fluctuations and the possibility of default [44]. Supposing that the covered event does not occur within the specified period of the earthquake bond, investors will be reimbursed their original amount and additional compensation for their exposure to catastrophic risk. In the event of an earthquake risk event with predetermined conditions, the SPV will compensate the sponsor, leading to partial or complete principal repayment to the investors [45,46]. The earthquake bond structure explained in this paragraph is visually summarized in Figure 2.

Determining trigger events is crucial in structuring earthquake bonds. There are three distinct categories of earthquake triggers: indemnity, index, and parametric [47]. An indemnity trigger pertains to the tangible financial damages experienced by the sponsor. Then, an industry index trigger uses the utilization of an index derived from estimations of losses provided by catastrophe-modeling firms, such as EQECAT, Applied Insurance Research (AIR) Worldwide, Property Claim Service (PCS), and Risk Management Solutions (RMS). Lastly, a parametric trigger relies on measurements of strength scale or depth of earthquake magnitude from designated data stations. In this paper, we model earthquake bond prices that utilize parametric triggers based on the strength scale.

3. Modeling Earthquake Bond Prices

3.1. Mathematical Notations and Assumptions

The following are the assumptions made in EB price modeling used in this research:

• (a). $\left[0, T\right]$ for a fixed $T>0$ is a real number interval representing continuous trading time;

• (b). $(\mathsf{\Omega },\mathcal{F},P)$ represents probability space, where $\mathsf{\Omega }$ is the set of states of the world, $\mathcal{F}$ is the sigma field of subsets of $\mathsf{\Omega }$, and $P$ is the probability measure on $\mathcal{F}$;

• (c). ${\mathcal{F}}_t\subset \mathcal{F}$ for $t\in \left[0,T\right]$ is increased filtration;

• (d). $\left\{N_t:t\in \left[0, T\right]\right\}$ is an inhomogeneous Poisson process with intensity $m_t={\int }_0^t{\lambda }_{s}ds>0$ representing the process of extreme earthquake frequency;

• (e). $\mu$ is the threshold value between extreme and nonextreme earthquakes;

• (f). $\left\{X_j:j=\mathrm{1,2},\dots ,N_t\right\}$ is a sequence of random variables with the domain $\left\{x\in \mathbb{R}:x\ge \mu \right\}$ representing the strength of extreme earthquakes. These are assumed to be independent and identically distributed (i.i.d.);

• (g). $\left\{M_t:t\in \left[0, T\right]\right\}$ represents the maximum strength process of an extreme earthquake;

• (h). $\left\{V_t:t\in \left[0, T\right]\right\}$ represents the EB price process involving region type, earthquake type, interest rate, extreme earthquake intensity, etc.;

• (i). $\left\{r_t:t\in \left[0,T\right]\right\}$ represents the force of interest rate process.

3.2. Force of Interest Rate Process

In general, the force of interest rates throughout the world has a nonnegative value. Only a few countries have negative values, e.g., Japan and Switzerland. Therefore, in this study, we assume that $\left\{r_t:t\in \left[0, T\right]\right\}$ is a nonnegative process and follows the mean-reverting square root process, which is formulated as follows [48]:

(1)$d{r}_t=\kappa \left(\theta -r_t\right)dt+{\sigma }_1\sqrt{r_t}d{W}_t,$

where $\kappa >0$ represents the mean reversion force measurement, $\theta >0$ represents the long-run mean of the force of interest, ${\sigma }_1>0$ represents the volatility parameter for the force of interest, and $\left\{W_t:t\in \left[0, T\right]\right\}$ represents the standard Wiener process according to physical measure $Q$ initiated at zero. In this process, the force of interest is guaranteed to be nonnegative because $2\kappa \theta \ge {\sigma }_1^2$ [49]. This process is also known as Cox–the Ingersoll–Ross model. If the model is applied to a country with a negative force of interest rates, several models can be used, e.g., the Vasicek and the Hull–White models.

The process in Equation (1) can be used for pricing zero-coupon bonds under the no-arbitrage assumption. In summary, the price of a zero-coupon bond with a maturity of $[t, T]$ is given as follows [34]:

(2)$E^Q\left(e^{-{\int }_t^T{r}_{s}ds}|{\mathcal{F}}_t\right)=B_{CIR}\left(t,T\right)=A\left(t,T\right)e^{-B\left(t,T\right)r_t},$

where

(3)$$\begin{gathered} A\left(t,T\right)={\left[\frac{2\gamma e^{\frac{\left(\kappa +\gamma \right)\left(T-t\right)}{2}}}{\left(\kappa +\gamma \right)\left(e^{\gamma \left(T-t\right)}-1\right)+2\gamma }\right]}^{\frac{2\kappa \theta }{{\sigma }_1^2}}, \\ B\left(t,T\right)=\frac{2\left(e^{\gamma \left(T-t\right)}-1\right)}{\left(\kappa +\gamma \right)\left(e^{\gamma \left(T-t\right)}-1\right)+2\gamma }, \\ \gamma =\sqrt{{\kappa }^2+2{\sigma }_1^2}. \end{gathered}$$

3.3. Extreme Earthquake Strength Process

The strength of extreme earthquakes in this study is represented as a sequence of random variables $\left\{X_j:j=\mathrm{1,2},\dots , N_t\right\}$ with domains $\left\{x\in \mathbb{R}:x\ge \mu \right\}$ that are i.i.d. The distribution function of $X_j$ is as follows:

(4)$F_X\left(x\right)=P\{X\le x|X>\mu \}=\frac{P\left\{\mu \le X\le x\right\}}{P\left\{X>\mu \right\}}.$

Equation (4) is difficult to determine because the “parent” distribution of $X_j$ is unknown. Therefore, we use the Pickands–Balkema–de Haan Theorem (See Salvadori et al. [50] and Tang and Yuan [13]) as an alternative to approach it. If there exist functions $𝒶\left(\mu \right)$ and $𝒷\left(\mu \right)$ with $𝒶\left(\mu \right)>0$, then when $\mu \to \infty$, $F\left[𝒶\left(\mu \right)\left(x-\mu \right)+𝒷\left(\mu \right)\right]$ converges to the generalized Pareto distribution as follows:

(5)$G_{\xi ,{\sigma }_2,\mu }\left(x\right)=\left\{\begin{array}{cc}1-{\left[1+\frac{\xi }{{\sigma }_2}\left(x-\mu \right)\right]}^{-\frac{1}{\xi }}& :\xi \ne 0\\ 1-e^{-\frac{x-\mu }{{\sigma }_2}}& :\xi =0\end{array}\right.,$

where $\xi \in \left(-\infty ,\infty \right)$ and ${\sigma }_2\in \left(0,\infty \right)$. When $\xi \ge 0$, the domain of $X_j$ is $\left\{x\in \mathbb{R}:x\ge \mu \right\}$, and when $\xi <0$, the domain of $X_j$ is $\left\{x\in \mathbb{R}:\mu \le x<\mu -\frac{{\sigma }_2}{\xi }\right\}$. In other words, this means that $F_X\left(x\right)$ with $x>\mu$ converges to $G_{\xi ,{\sigma }_2,\mu }\left(x\right)$ [51,52], or

(6)$F_X\left(x\right)=P\{X\le x|X>\mu \}=\frac{P\left\{\mu \le X\le x\right\}}{P\left\{X>\mu \right\}}\approx \left\{\begin{array}{cc}1-{\left[1+\frac{\xi }{{\sigma }_2}\left(x-\mu \right)\right]}^{-\frac{1}{\xi }}& :\xi \ne 0\\ 1-e^{-\frac{x-\mu }{{\sigma }_2}}& :\xi =0\end{array}\right..$

3.4. Maximum Strength of Extreme Earthquakes Process

Earthquakes occur all the time, so the event frequency is challenging to calculate. In other words, the intensity of the event in each time interval is infinite. Therefore, as an alternative, measurements of the risk of earthquake strength focus on extreme types of earthquakes only because their frequency can be calculated. Extreme earthquakes are defined as earthquakes with a strength greater than a certain threshold value.

Let $\mu$ be the threshold value between extreme and inextreme earthquake strengths. The maximum strength of an extreme earthquake is designed into a compound distribution consisting of the following:

• (a). Frequency of extreme earthquake events.

This is represented as an inhomogeneous Poisson process $\left\{N_t:t\in \left[0,T\right]\right\}$ with intensity $m_t={\int }_0^t{\lambda }_{s}ds>0$ and a probability function of [53,54]

(7)$\mathrm{P}\left\{N_t=n\right\}=\frac{m_t^n}{n!}{e}^{-m_t}:n=0,1,2,\dots ;$

• (b). The strength of an extreme earthquake.

This is represented as a sequence of random variables $\left\{X_j:j=\mathrm{1,2},\dots ,N_t\right\}$ with domain $\left\{x\in \mathbb{R}:x\ge \mu \right\}$ whose i.i.d. $G_{\xi ,{\sigma }_2,\mu }$.

The maximum value of extreme earthquake strength up to time $t$, denoted by $M_t$, is modeled as

(8)$M_t=\max \left\{X_j:j=\mathrm{1,2},\dots ,N_t\right\}.$

In more detail, $\left\{X_j:j=\mathrm{1,2},\dots ,N_t\right\}$ and $\left\{N_t:t\in \left[0,T\right]\right\}$ are assumed to be independent. In other words, the frequency of extreme earthquakes does not affect their strength, and vice versa. The distribution function of $M_t$ is

(9)$F_{M_t}\left(x\right)=\left\{\begin{array}{cc}{e}^{-m_t{\left(1+\frac{\xi }{{\sigma }_2}\left(x-\mu \right)\right)}^{-\frac{1}{\xi }}}& :\xi \ne 0\\ e^{-m_t{e}^{\left(-\frac{x-\mu }{{\sigma }_2}\right)}}& :\xi =0\end{array}\right..$

Proof. 

See Appendix A. □

3.5. Earthquake Bond Pricing Models

We use a risk-neutral pricing measure to model EB prices. We follow Cox and Pedersen [55], whereby the original distribution characteristics of the maximum strength of the extreme earthquake process are assumed to be retained after the conversion of the physical probability measure $P$ to the neutral risk measure $Q$. Under the neutral risk pricing measure $Q$, events that depend on financial variables and events that depend on earthquake risk variables do not influence each other (for more details, see Cox and Pedersen [55]).

Suppose a bond is purchased at time $t$ and matures at time $T$. In a zero-coupon EB, the investor will receive a payoff of $R_T$. This $R_T$ value is the payoff of the EB redemption value, which depends on the maximum strength of the extreme earthquake up to time $T$. In this study, $R_T$ is used as a piecewise function, which is mathematically expressed as follows:

(10)$R_T=\left\{\begin{array}{cc}{\eta }_{1}K& :0\le M_{[t,T]}<d_1\hfill \\ {\eta }_{2}K& :d_1\le M_{[t,T]}<d_2\hfill \\ ⋮& :\hspace{1em}\hspace{1em}\hspace{1em} ⋮\hfill \\ {\eta }_{\mathcal{S}-1}K& :d_{\mathcal{S}-2}\le M_{[t,T]}<d_{\mathcal{S}-1}\hfill \\ {\eta }_{\mathcal{S}}K& :M_{[t,T]}\ge d_{\mathcal{S}-1}\hfill \end{array}\right.,$

where $K>0$ represents the redemption value of EB, $\left\{{\eta }_1,{\eta }_2,\dots ,{\eta }_S:0\le {\eta }_S<{\eta }_{S-1}<\cdots <{\eta }_1\le 1,{\eta }_1+{\eta }_2+\cdots +{\eta }_S=1\right\}$ is a fixed value that represents the proportion of payoff $K$, $\left\{d_1,d_2,\dots ,d_{\mathcal{S}-1}\right\}$ is an increasing sequence of positive real numbers that represent the strength of the earthquake, and $M_{\left[t, T\right]}$ is the maximum strength of the extreme earthquake from time $t$ to time $T$. Using expectations under the risk-neutralized pricing measure $Q$, the zero-coupon price EB is expressed as follows:

(11)$V_t=E^Q\left(e^{-{\int }_t^T{r}_{s}ds}{R}_T|{\mathcal{F}}_t\right)=B_{CIR}\left(t,T\right)K\left[{\eta }_{\mathcal{S}}+\sum _{𝓂=1}^{\mathcal{S}-1}({\eta }_{𝓂}-{\eta }_{𝓂+1})P(M_{\left[t,T\right]}\le {𝒹}_{𝓂})\right],$

where

(12)$P\left(M_{\left[t,T\right]}\le {𝒹}_{𝓂}\right)=\left\{\begin{array}{cc}{e}^{-m_{\left[t,T\right]}{\left(1+\frac{\xi }{{\sigma }_2}\left({𝒹}_{𝓂}-\mu \right)\right)}^{-\frac{1}{\xi }}}& :\xi \ne 0\\ e^{-m_{\left[t,T\right]}{e}^{\left(-\frac{{𝒹}_{𝓂}-\mu }{{\sigma }_2}\right)}}& :\xi =0\end{array}\right.,$

and

(13)$$\begin{gathered} B_{CIR}\left(t,T\right)=A\left(t,T\right)e^{-B\left(t,T\right)r_t}, \\ A\left(t,T\right)={\left[\frac{2\gamma e^{\frac{\left(\kappa +\gamma \right)\left(T-t\right)}{2}}}{\left(\kappa +\gamma \right)\left(e^{\gamma \left(T-t\right)}-1\right)+2\gamma }\right]}^{\frac{2\kappa \theta }{{\sigma }_1^2}}, \\ B\left(t,T\right)=\frac{2\left(e^{\gamma \left(T-t\right)}-1\right)}{\left(\kappa +\gamma \right)\left(e^{\gamma \left(T-t\right)}-1\right)+2\gamma }, \\ \gamma =\sqrt{{\kappa }^2+2{\sigma }_1^2}. \end{gathered}$$

Proof. 

See Appendix B. □

Next is coupon-paying EB (CPEB) price modeling. At time $T$, investors in this type of EB will receive a payoff of $R_T$. The amount of $R_T$ consists of the redemption value and coupon, both of which depend on the maximum strength of the extreme earthquake up to time $T$. $R_T$ is reconstructed as a piecewise function, which is mathematically expressed as follows:

(14)$R_T=\left\{\begin{array}{cc}{\eta }_1(K+C)& :0\le M_{\left[t,T\right]}<{𝒹}_1 \\ {\eta }_2(K+C)& :{𝒹}_1\le M_{\left[t,T\right]}<{𝒹}_2 \\ ⋮& : ⋮ \\ {\eta }_{\mathcal{S}-1}(K+C)& :{𝒹}_{\mathcal{S}-2}\le M_{\left[t,T\right]}<{𝒹}_{\mathcal{S}-1}\\ {\eta }_{\mathcal{S}}(K+C)& :M_{\left[t,T\right]}\ge d_{\mathcal{S}-1} \end{array},\right.$

where $C$ represents the coupon. Using expectation under the risk-neutralized pricing measure $Q$, the CPEB price model is expressed as follows:

(15)$V_t=E^Q\left(e^{-{\int }_t^T{r}_{s}ds}{R}_T|{\mathcal{F}}_t\right)=B_{CIR}\left(t,T\right)\left(K+C\right)\left[{\eta }_{\mathcal{S}}+\sum _{𝓂=1}^{\mathcal{S}-1}\left({\eta }_{𝓂}-{\eta }_{𝓂+1}\right)P\left(M_{\left[t,T\right]}\le {𝒹}_{𝓂}\right)\right],$

where

(16)$P\left(M_{\left[t,T\right]}\le {𝒹}_{𝓂}\right)=\left\{\begin{array}{cc}{e}^{-m_{\left[t,T\right]}{\left(1+\frac{\xi }{{\sigma }_2}\left({𝒹}_{𝓂}-\mu \right)\right)}^{-\frac{1}{\xi }}}& :\xi \ne 0\\ e^{-m_{\left[t,T\right]}{e}^{\left(-\frac{{𝒹}_{𝓂}-\mu }{{\sigma }_2}\right)}}& :\xi =0\end{array}\right.,$

and

(17)$$\begin{gathered} B_{CIR}\left(t,T\right)=A\left(t,T\right)e^{-B\left(t,T\right)r_t}, \\ A\left(t,T\right)={\left[\frac{2\gamma e^{\frac{\left(\kappa +\gamma \right)\left(T-t\right)}{2}}}{\left(\kappa +\gamma \right)\left(e^{\gamma \left(T-t\right)}-1\right)+2\gamma }\right]}^{\frac{2\kappa \theta }{{\sigma }_1^2}}, \\ B\left(t,T\right)=\frac{2\left(e^{\gamma \left(T-t\right)}-1\right)}{\left(\kappa +\gamma \right)\left(e^{\gamma \left(T-t\right)}-1\right)+2\gamma }, \\ \gamma =\sqrt{{\kappa }^2+2{\sigma }_1^2}. \end{gathered}$$

Proof. 

See Appendix C. □

4. Model Experiment to Actual Data

4.1. Data Description

Indonesia has the second highest frequency of earthquakes in the world after Japan. This happens because the region is geologically located above four active plate lines: Eurasia, Indo-Australia, the Philippines, and the Pacific. However, under these conditions, Indonesia has never issued EBs [56,57,58]. Therefore, the model in this research is applied to earthquake data in Indonesia.

We obtained data on the strength and frequency of extreme earthquakes from the National Disaster Management Authority (BNPB) of the Republic of Indonesia from 2008 to 2021. We use a threshold value between extreme and inextreme earthquake strengths of $\mu =5$ because earthquakes with strengths above this value cause significant damage. Data on the strength of extreme earthquakes were collected daily from 1 November 2008 to 31 December 2021 (4808 days). Meanwhile, data on the frequency of extreme earthquakes were collected annually from 2008 to 2021. Visually, data on extreme earthquake strength and frequency are given in Figure 3 and Figure 4, respectively. Statistical descriptions of both are given in Table 2. Table 2 shows that the means of strength and frequency of extreme earthquakes in Indonesia are 5.3452 $M_L$ (Richter scale) and 258.5 extreme earthquakes per year, respectively. The standard deviations are 0.3794 $M_L$ and 119.16488 extreme earthquakes per year, respectively. In one year, there are 596 extreme earthquakes on average, while the highest extreme earthquake strength is 7.9 $M_L$.

4.2. Extreme Earthquake Strength Distribution

In this research, the theoretical distribution that will be matched to the distribution of extreme earthquake strength data is the GPD given by Equation (5). There are two types of GPD used in this equation: when $\xi \ne 0$, which we call GPD Type I, and when $\xi =0$, which we call GPD Type II. Parameter estimation for the two types of GPD was conducted first using the MLE method. The results of the parameter estimate from the two distributions are given in Table 3.

Then, the next step is checking the fit between the two types of GPD and the data distribution. This was conducted visually and theoretically. Visually, it was conducted using a probability–probability plot (PP-Plot). The data distribution follows the GPD visually if both scatters on the Cartesian diagram are around the line with one gradient [59]. The PP-Plot of the data distribution and the two types of GPD are given in Figure 5.

Figure 5 shows that the scatter between the data distribution with GPD Type I is closer around the one gradient line colored red than the scatter between the data distribution with GPD Type II. This indicates that the data distribution follows GPD Type I more closely than GPD Type II. Next is a theoretical examination using the Anderson–Darling, Kolmogorov–Smirnov, and Chi-Squared tests, with a significant value of 0.05 and a sample size 3468. For these three goodness-of-fit tests, if the test value is smaller than the critical value, then the data distribution follows the specified theoretical distribution, and vice versa [60,61,62]. Briefly, the results of the fit check are given in Table 4.

Table 4 shows that the GPD Type I test values in each statistical test are smaller than the critical values, while all GPD Type II test values in each statistical test are not smaller than the critical values. This indicates that the data distribution follows GPD Type I. Therefore, GPD Type I was chosen to describe the data distribution.

4.3. Extreme Earthquake Frequency Process

In this study, the process of determining the annual frequency of extreme earthquakes is modeled using the inhomogeneous Poison process $\left\{N_t:t\in \left[0, T\right]\right\}$ with intensity $m_t={\int }_0^t{\lambda }_{s}ds>0$. In Figure 4, the data pattern of annual extreme earthquake intensity appears to tend towards being cyclical and periodic. Therefore, we use a model that can accommodate this, namely the Vere–Jones model [63]. The model is stated as follows:

(18)${\lambda }_t=a+b\sin (t+c)+{\lambda }_c,$

where $a$, $b\in \mathbb{R}$. Equation (18) is hereafter called model A. Model A can be discretized. The method is to discretize the interval $\left[0, T\right]$ into $\left\{{𝓉}_0=0,{𝓉}_1,\dots ,{𝓉}_z=T\right\}$, where $z$ represents the number of steps. The discretization results of Model A are expressed as follows:

${\lambda }_{{𝓉}_{k+1}}=a+b\sin \left({𝓉}_{k+1}+{𝓉}_k\right)+{\lambda }_{{𝓉}_k},$

(19)$\mathsf{\Delta }{\lambda }_{{𝓉}_{k+1}}=a+b\sin \left({𝓉}_{k+1}+{𝓉}_k\right),$

where $\mathsf{\Delta }{\lambda }_{{𝓉}_{k+1}}={\lambda }_{{𝓉}_{k+1}}-{\lambda }_{{𝓉}_k}$ with $k=1, 2,\dots , z-1$. The MLE method is used to estimate parameters $a$ and $b$. Briefly, the estimation results $\widehat{a}=-8.13277$ and $\widehat{b}=21.37195$. In addition, we also compare the model in Equation (19) with the autoregressive integrated moving average (ARIMA) (1, 1, 1) model, as follows [64]:

(20)$\mathsf{\Delta }{\lambda }_{{𝓉}_{k+1}}={\alpha }_1\mathsf{\Delta }{\lambda }_{{𝓉}_k}+{\beta }_1{ϵ}_{{𝓉}_k}+{ϵ}_{{𝓉}_{k+1}},$

where ${\alpha }_1$ $\in$ $\mathbb{R}$ represents the autoregressive coefficient, ${\beta }_1$ $\in$ $\mathbb{R}$ represents the moving average coefficient, and ${ϵ}_{{𝓉}_k}~i.i.d. N(0,{\sigma }^2)$. Equation (20) is hereafter called model B. To estimate the parameters ${\alpha }_1$ and ${\beta }_1$, the MLE method is used. Briefly, the estimation results ${\widehat{\alpha }}_1=-0.8610$ and ${\widehat{\beta }}_1=0.1734$. We also considered using a constant intensity $\widehat{\lambda }=258.5$. This is referred to as model C. For model selection, the accuracy measure used is a mean absolute percentage error (MAPE). The results of these accuracy measurements produce a MAPE of 57 per cent for model A, 30 per cent for model B, and 37 per cent for model C. Based on the criteria of time-series model accuracy from Lewis [65], model A is said to be inaccurate, while models B and C are called reasonable. Because the MAPE of model B is smaller than that of model C, the model used to estimate the annual intensity of extreme earthquakes in this experiment is model B, ARIMA (1, 1, 1).

4.4. Force of Interest Rate Process

The data used to apply the CIR model in this research are annual interest rate data from Indonesia from 1990 to 2021. The data are converted into the force of interest with the formula $r_t=\ln \left(1+{𝒾}_t\right)$, where ${𝒾}_t$ is the interest rate at time $t$. Data were obtained from Bank Indonesia, the Central Bank in the country. First, we discretized Equation (1) to make it easier to estimate the parameters. We discretized the interval $\left[0, T\right]$ into $\left\{t_0=0,t_1,\dots ,t_{𝓃}=T\right\}$, where $𝓃$ represents the number of steps. The discretization result of Equation (1) is expressed as follows [66]:

$\mathsf{\Delta }{r}_{t_{i+1}}=\kappa \left(\theta -r_{t_i}\right)\mathsf{\Delta }{t}_{i+1}+{\sigma }_1\sqrt{r_{t_i}\mathsf{\Delta }{t}_{i+1}}{\epsilon }_{t_i},$

(21)$\frac{\mathsf{\Delta }{r}_{t_{i+1}}}{\sqrt{r_{t_i}}}=\kappa \theta \frac{\mathsf{\Delta }{t}_{i+1}}{\sqrt{r_{t_i}}}-\kappa \sqrt{r_{t_i}}\mathsf{\Delta }{t}_{i+1}+{\sigma }_1\sqrt{\mathsf{\Delta }{t}_{i+1}}{\epsilon }_i,$

where $i=1,2,\dots ,n-1$, $\Delta r_{t_{i+1}}=r_{t_{i+1}}-r_{t_i}$, $\Delta t_{i+1}=t_{i+1}-t_i$, and ${\epsilon }_i$ represents i.i.d. $N(0,\mathsf{\Delta }{t}_{i+1})$. Parameter estimation of Equation (21) was conducted 15,000 times using the MLE method. Briefly, the results of these parameter estimates are $\widehat{\kappa }=0.493096$, $\widehat{\theta }=0.0255701$, and ${\widehat{\sigma }}_1=0.002278$. We also compare the CIR model with the constant interest rate model $r=0.098015$. We call this the constant interest rate model. Next, the accuracy of both models was measured using MAPE. In short, the MAPE obtained was 31 per cent for the CIR model and 54 per cent for the constant interest rate model. According to Lewis [65], the accuracy of the CIR model was reasonable, while the accuracy of the constant interest rate model was low. Therefore, the CIR model was chosen to estimate the annual interest rate in the model experiments in this study.

4.5. Estimating Earthquake Bond Prices

To determine the estimated price of zero-coupon EBs (ZCEBs) and coupon-paying EBs (CPEBs), the values of some variables must be determined first. The values are given in Table 5.

Table 5 shows that an EB has a redemption value of $K=1$ IDR, a coupon of $C=0.1$ IDR, and a term of 2 years, namely, from 1 January 2022 to 31 December 2023. By using the ARIMA (1, 1, 1) model, we iteratively obtained that the intensities of extreme earthquakes in 2022 and 2023 were ${\widehat{\lambda }}_{15}=261.2826$ and ${\widehat{\lambda }}_{16}=264.5583$ extreme earthquakes per year, respectively. Hence, $m_{\left[t,T\right]}={\int }_t^T{\lambda }_{s}ds\approx \sum _{𝓀=15}^{16}{\widehat{\lambda }}_{𝓀}=525.8409$. Then, using the force of interest in 2021, $r_t\approx r_{34}=0.0344014$, we obtained $B_{CIR}\left(t, T\right)\approx 0.98112$. Then, the payoffs of ZCEB and CPEB at time $T$ were as follows:

(22)$R_T=\left\{\begin{array}{cc}1& :0\le M_{\left[t,T\right]}<5\hfill \\ 0.875& :5\le M_{\left[t,T\right]}<6\hfill \\ 0.75& :6\le M_{\left[t,T\right]}<7\hfill \\ 0.625& :7\le M_{\left[t,T\right]}<8\hfill \\ 0.5& :M_{\left[t,T\right]}\ge 8\hfill \end{array}\right.$

and

(23)$R_T=\left\{\begin{array}{cc}1.1& :0\le M_{\left[t,T\right]}<5\hfill \\ 0.875\left(1.1\right)& :5\le M_{\left[t,T\right]}<6\hfill \\ 0.75\left(1.1\right)& :6\le M_{\left[t,T\right]}<7\hfill \\ 0.625\left(1.1\right)& :7\le M_{\left[t,T\right]}<8\hfill \\ 0.5\left(1.1\right)& :M_{\left[t,T\right]}\ge 8\hfill \end{array}\right.,$

respectively. The selection of the interval $M_{\left[t, T\right]}$ at each step is based on the classification of earthquakes in Indonesia based on their strength, as follows:

• (a). $0\le M_{\left[t, T\right]}<5$ represents the classification of earthquakes with micro to light strength;

• (b). $5\le M_{\left[t, T\right]}<6$ represents the classification of an earthquake with moderate strength;

• (c). $6\le M_{\left[t, T\right]}<7$ represents the classification of an earthquake with strong strength;

• (d). $7\le M_{\left[t,T\right]}<8$ represents the earthquake classification with major strength;

• (e). $M_{\left[t, T\right]}\ge 8$ represents the classification of an earthquake with great strength.

Then, the set of EB payoff proportions can be adjusted to suit investor tolerance. In this study, investors only want to lose half of their payoffs if the maximum strength of the earthquake is greater than or equal to 8 $M_L$ and receive the whole amount if the maximum strength of the earthquake is less than 5 $M_L$. Meanwhile, for other earthquake strength intervals, the proportion of payoffs to EB is determined by an arithmetic series with a difference of 0.125. With the specified variable values, the estimated prices of ZCEB and CPEB in Indonesia are IDR 0.5595 and IDR 0.6155, respectively.

5. Discussion

Returns from ZCEB and CPEB are as follows:

(24)$R_{ZCEB}=\left\{\begin{array}{cc}\frac{1-0.5595}{1}& :0\le M_{\left[t,T\right]}<5\hfill \\ \frac{0.875-0.5595}{0.875}& :5\le M_{\left[t,T\right]}<6\hfill \\ \frac{0.75-0.5595}{0.75}& :6\le M_{\left[t,T\right]}<7\hfill \\ \frac{0.625-0.5595}{0.625}& :7\le M_{\left[t,T\right]}<8\hfill \\ \frac{0.5-0.5595}{0.5}& :M_{\left[t,T\right]}\ge 8\hfill \end{array}\right.=\left\{\begin{array}{cc}0.4405& :0\le M_{\left[t,T\right]}<5\hfill \\ 0.3605& :5\le M_{\left[t,T\right]}<6\hfill \\ 0.2539& :6\le M_{\left[t,T\right]}<7\hfill \\ 0.1047& :7\le M_{\left[t,T\right]}<8\hfill \\ -0.1191& :M_{\left[t,T\right]}\ge 8\hfill \end{array}\right.$

and

(25)$R_{CPEB}=\left\{\begin{array}{cc}\frac{1.1-0.6155}{1.1}& :0\le M_{\left[t,T\right]}<5\hfill \\ \frac{0.875\left(1.1\right)-0.6155}{0.875\left(1.1\right)}& :5\le M_{\left[t,T\right]}<6\hfill \\ \frac{0.75\left(1.1\right)-0.6155}{0.75\left(1.1\right)}& :6\le M_{\left[t,T\right]}<7\hfill \\ \frac{0.625\left(1.1\right)-0.6155}{0.625\left(1.1\right)}& :7\le M_{\left[t,T\right]}<8\hfill \\ \frac{0.5\left(1.1\right)-0.6155}{0.5\left(1.1\right)}& :M_{\left[t,T\right]}\ge 8\hfill \end{array}\right.=\left\{\begin{array}{cc}0.4405& :0\le M_{\left[t,T\right]}<5\hfill \\ 0.3605& :5\le M_{\left[t,T\right]}<6\hfill \\ 0.2539& :6\le M_{\left[t,T\right]}<7\hfill \\ 0.1047& :7\le M_{\left[t,T\right]}<8\hfill \\ -0.1191& :M_{\left[t,T\right]}\ge 8\hfill \end{array}\right.,$

respectively. Equations (24) and (25) show that with the same payoff proportions $\left\{{\eta }_1,{\eta }_2,\dots ,{\eta }_5\right\}$ and intervals $M_{\left[t, T\right]}$, the returns from ZCEB and CPEB appear to be the same for each possibility. Thus, the means of returns of both are also the same. The maximum return from EBs based on Equations (24) and (25) is 44 per cent. This amount is very large compared to traditional bond securities. Then, the smallest return from EBs was deemed to be $-11.8$ per cent. The means of returns from ZCEB and CPEB received by investors are as follows:

(26)$$\begin{gathered} E^P\left(R_{ZCEB}|{\mathcal{F}}_t\right)=E^P\left(R_{CPEB}|{\mathcal{F}}_t\right), \\ =E^P(0.440791_{0\le M_{\left[t,T\right]}<5}+0.360901_{5\le M_{\left[t,T\right]}<6}+\dots -0.118421_{M_{\left[t,T\right]}\ge 8}|{\mathcal{F}}_t), \\ =0.44079P\left(0\le M_{\left[t,T\right]}<5\right)+0.36090P\left(5\le M_{\left[t,T\right]}<6\right)+\dots -0.11842P\left(M_{\left[t,T\right]}\ge 8\right), \\ =0.00494. \end{gathered}$$

Next, we examine the effect of the EB term on its price. Rationally, the EB term has a negative relationship with its price, whereby the longer the EB term, the lower the price, and vice versa. We will validate this through our model. Visually, the relationship between the EB term and its price is given in Figure 6. This shows that the EB term is inversely proportional to its price. Therefore, the model can describe this according to reality.

The effect of the stochastic force of interest on EB prices is analyzed next by examining the difference between EB prices with the stochastic and constant force of interest assumptions for a term of one to five years. The constant interest rate used is 6 per cent according to actual conditions in 2021. This difference in EB prices with two assumptions is visually shown in Figure 7. It shows significant differences in EB prices between the two assumptions, with ZCEB and CPEB prices ranging between 3.12 and 9.21 per cent and 3.43 and 10.13 per cent, respectively.

The analysis in Section 4.2 reveals that GPD Type I is more fit than GPD Type II for describing the distribution of extreme earthquake strength data in Indonesia. The study also examines the effect of the fit on EB prices by examining the difference between EB prices using two GPD types for a term of one to five years. The difference in EB prices is visually shown in Figure 8. It shows significant differences in EB prices, with ZCEB and CPEB prices ranging between $-7.01$ and $-6.28$ per cent and $-7.71$ and $-7.02$ per cent, respectively. Therefore, the careful selection of the most appropriate type of GPD is crucial.

Based on the analysis in Section 4.3, inconstant extreme earthquake intensity models are more accurate than constant ones in describing annual extreme earthquake intensity data. The effect of this accuracy on EB prices is also analyzed by examining the difference between EB prices using these inconstant and constant assumptions for a term of one to five years. The difference in EB prices is visually shown in Figure 9. Figure 9 shows significant differences in EB prices, with ZCEB and CPEB prices ranging between $-0.09$ and $-0.06$ per cent and $-0.10$ and $-0.07$ per cent, respectively. Therefore, it is crucial make the choice using assumptions to yield a more accurate annual intensity forecast.

6. Conclusions

This research presents an EB pricing model involving inconstant event intensity and the maximum strength of extreme earthquakes under the risk-neutral pricing measure. Focusing on extreme earthquakes simplifies the modeling process and data collection and computing time compared to considering the infinite frequency of earthquakes occurring over a continuous time interval. The inconstant intensity of the event is accommodated by an inhomogeneous Poisson process, while the maximum strength is modeled using extreme value theory (EVT). Then, the model is applied to earthquake data in Indonesia, the country with the second highest frequency of earthquakes worldwide. Finally, the variable sensitivities of EB prices are also analyzed.

The results of the sensitivity analysis of the variables show that EB can provide a positive return of almost 50 per cent. Meanwhile, it is also possible to obtain the largest negative return, but the quantity is not as much as the maximum positive return. Then, the inconstant extreme earthquake intensity and force of interest rate have crucial implications for EB prices. There is a significant difference between EB prices with the assumption of inconstancy and constancy in two variables. This indicates that it is important to use the inconstant extreme earthquake intensity and force of interest rate. Then, the selection of GPD Types I and II that fit the extreme earthquake strength data distribution also has an important effect on the EB price, and there is also a significant difference in EB prices based on the fits of the two GPD types. Therefore, the careful selection of the most appropriate type of GPD is crucial.

The maximum strength risk model for extreme earthquakes can be utilized by practitioners, disaster management agencies, or geological agencies to measure the risk of such disasters in a region, and the EB price models can be used by the issuer in setting EB prices. Then, the conducted experiments can guide how to use the models designed in this research. Investors can use the sensitivity of the variables to EB prices to choose EBs that suit their risk tolerance. As a suggestion for future research, the model in this study can be developed to include multiple triggers. An additional trigger that can be used is the indemnity index. The involvement of this index ensures that the severity of the earthquake can also be measured financially. Then, the country’s territory can be further decomposed into smaller administrative regions. Accommodating differences in earthquake risk characteristics through this decomposition can help make earthquake bond prices more reasonable.

Author Contributions

Conceptualization, R.A.I., S., H.N. and R.I.I.; methodology, R.A.I.; software, R.A.I.; validation, S., H.N. and R.A.I.; formal analysis, R.A.I.; investigation, R.A.I.; resources, S.; data curation, R.A.I.; writing—original draft preparation, R.A.I.; writing—review and editing, S., H.N. and R.I.I.; visualization, R.A.I.; supervision, S., H.N. and R.I.I.; project administration, S.; funding acquisition, S. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data are contained within the article.

Acknowledgments

Thanks are conveyed to Universitas Padjadjaran, who provided the Padjadjaran Doctoral Program Scholarship.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Average of loss from several types of extreme disasters worldwide.

Figure 2. Simple framework of earthquake bond.

Figure 3. Extreme earthquake strength from 1 November 2008 to 31 December 2021 in Indonesia.

Figure 4. Annual frequency of extreme earthquakes from 2008 to 2021 in Indonesia.

Figure 5. Scatter between data distribution with GPD Type I (a) and Type II (b).

Figure 6. Effect of EB term on its price.

Figure 7. The difference between EB prices with the assumption of stochastic and constant force of interest rates.

Figure 8. The difference between EB prices based on the use of GPD Types I and II.

Figure 9. The difference between EB prices using inconstant and constant extreme earthquake intensities.

Table 1

History of earthquake bond issuance until 2022.

Month of Issues	Country	Sponsor	Special Purpose Vehicles	Term (Years)	Contingency Fund(USD Million)
May 2006	Mexico	FONDEN	CAT-Mex Ltd.	3	160
October 2009	Mexico	FONDEN	MultiCAT Mexico Ltd.	2	290
August 2017	Mexico	FONDEN	IBRD CAR 113	2	120
February 2018	Chile	Chile	IBRD CAR 116	3	500
February 2018	Colombia	Colombia	IBRD CAR 117	3	400
February 2018	Mexico	FONDEN	IBRD CAR 118-119	2	260
February 2018	Peru	Peru	IBRD CAR 120	2	200
November 2019	Philippines	Philippines	IBRD CAR 123	3	75

Table 2

Statistical description of the strength and frequency of extreme earthquakes in Indonesia.

Statistical Description	Value
Mean of Strength ($M_L$)	5.3452
Maximum of Strength ($M_L$)	7.9
Standard Deviation of Strength ($M_L$)	0.3794
Mean of Frequency (Earthquake per year)	258.5
Maximum of Frequency (Earthquake per year)	596
Standard Deviation of Frequency (Earthquake per year)	119.16488

Table 3

GPD parameter estimation results.

GPD Type	Parameter Estimator(s)
I	$\widehat{\xi }=0.100956$ and ${\widehat{\sigma }}_2=0.3106285$
II	${\widehat{\sigma }}_2=0.3452422$

Table 4

Distribution fit test results.

Statistical Test	Critical Value (for 0.05 Significant Value)	Test Value
		GPD Type I	GPD Type II
Anderson–Darling	2.5018	1.3707	7.8682
Kolmogorov–Smirnov	0.0231	0.0155	0.0389
Chi-Squared	19.6750	15.8630	47.6670

Table 5

Variable values in EB price estimation.

Variable	Value
$t$	1 January 2022
$T$	31 December 2023
$K$	1 IDR
$C$	0.1 IDR

Appendix A

The event $\{M_t\le x\}$ occurs if $n$ extreme earthquakes occur, $n=0, 1, 2,\dots ,$ and the maximum strength of the $n$ extreme earthquakes is less than or equal to $x$. Therefore, the event $\{M_t\le x\}$ can be represented as a union of the mutually exclusive events $\{N_t=n,M_t\le x\}$, or $$\left\{M_t\le x\right\}=\bigcup _{n=0}^{\infty }\{N_t=n,M_t\le x\}.$$ Thus, $$\begin{array}{l}{F}_{M_t}\left(x\right)=P\left\{M_t\le x\right\},\\ ={\displaystyle \sum _{n=0}^{\infty }}P\{N_t=n,M_t\le x\},\\ ={\displaystyle \sum _{n=0}^{\infty }}P\{N_t=n\}P\left\{M_t\le x|N_t=n\right\}.\end{array}$$ Then, the event $\left\{M_t\le x\right\}$ for a fixed $n$ is equivalent to the event $\left\{X_1\le x,\dots ,X_n\le x|X_1>\mu ,\dots ,X_n>\mu \right\}$ such that $$\begin{array}{l}{F}_{M_t}\left(x\right)={\displaystyle \sum _{n=0}^{\infty }}P\{N_t=n\}P\left\{X_1\le x,\dots ,X_n\le x|X_1>\mu ,\dots ,X_n>\mu \right\},\\ ={\displaystyle \sum _{n=0}^{\infty }}P\left\{N_t=n\right\}\frac{P\left\{X_1\le x,\dots ,X_n\le x,X_1>\mu ,\dots ,X_n>\mu \right\}}{P\left\{X_1>\mu ,\dots ,X_n>\mu \right\}},\\ ={\displaystyle \sum _{n=0}^{\infty }}P\left\{N_t=n\right\}\frac{P\left\{\mu \le X_1\le x,\dots ,\mu <X_n\le x\right\}}{P\{X_1>\mu ,\dots ,X_n>\mu \}},\\ ={\displaystyle \sum _{n=0}^{\infty }}P\left\{N_t=n\right\}\frac{P\left\{\mu \le X_1\le x\right\}\dots P\{\mu <X_n\le x\}}{P\left\{X_1>\mu \right\}\dots P\left\{X_n>\mu \right\}},\\ ={\displaystyle \sum _{n=0}^{\infty }}P\left\{N_t=n\right\}{\left(\frac{P\left\{\mu \le X\le x\right\}}{P\left\{X>\mu \right\}}\right)}^n,\\ ={\displaystyle \sum _{n=0}^{\infty }}\frac{m_t^n}{n!}{e}^{-m_t}{\left(G_{\xi ,{\sigma }_2,\mu }\left(x\right)\right)}^n,\\ ={\displaystyle \sum _{n=0}^{\infty }}\frac{{\left(m_t{G}_{\xi ,{\sigma }_2,\mu }\left(x\right)\right)}^n}{n!}{e}^{-m_t},\\ ={\displaystyle \sum _{n=0}^{\infty }}\frac{{\left(m_t{G}_{\xi ,{\sigma }_2,\mu }\left(x\right)\right)}^n}{n!}{e}^{-m_t}\frac{e^{-m_t\left(G_{\xi ,{\sigma }_2,\mu }\left(x\right)-1\right)}}{e^{-m_t\left(G_{\xi ,{\sigma }_2,\mu }\left(x\right)-1\right)}},\\ =\frac{1}{e^{-m_t\left(G_{\xi ,{\sigma }_2,\mu }\left(x\right)-1\right)}}{\displaystyle \sum _{n=0}^{\infty }}\frac{{\left(m_t{G}_{\xi ,{\sigma }_2,\mu }\left(x\right)\right)}^n}{n!}{e}^{-m_t{G}_{\xi ,{\sigma }_2,\mu }\left(x\right)},\\ =e^{m_t\left(G_{\xi ,{\sigma }_2,\mu }\left(x\right)-1\right)},\\ =\left\{\begin{array}{cc}{e}^{m_t\left[1{-\left(1+\frac{\xi }{{\sigma }_2}\left(x-\mu \right)\right)}^{-\frac{1}{\xi }}-1\right]}& :\xi \ne 0\\ e^{m_t\left[1-e^{\left(-\frac{x-\mu }{{\sigma }_2}\right)}-1\right]}& :\xi =0\end{array}\right.,\\ =\left\{\begin{array}{cc}{e}^{-m_t{\left(1+\frac{\xi }{{\sigma }_2}\left(x-\mu \right)\right)}^{-\frac{1}{\xi }}}& :\xi \ne 0\\ e^{-m_t{e}^{\left(-\frac{x-\mu }{{\sigma }_2}\right)}}& :\xi =0\end{array}\right..\end{array}$$

Appendix B

Let $\left[0,T\right]$ for finite positive real number $T$ be a continuous trading time. Let financial market and earthquake risk variables be modeled within filtered probability triples $\left({\mathsf{\Omega }}^{\left(1\right)}, {\mathcal{F}}^{\left(1\right)},P^{\left(1\right)}\right)$ and $({\mathsf{\Omega }}^{\left(2\right)}, {\mathcal{F}}^{\left(2\right)},P^{\left(2\right)})$, respectively. Then, let $\left(\mathsf{\Omega },\mathcal{F},P\right)$ be the probability triple for the full model, where:

• $\mathsf{\Omega }={\mathsf{\Omega }}^{\left(1\right)}\times {\mathsf{\Omega }}^{\left(2\right)}$ represents the sample space for the full model;

• $\mathcal{F}={\mathcal{F}}^{\left(1\right)}\times {\mathcal{F}}^{\left(2\right)}$ represents the sigma field of subsets of $\mathsf{\Omega }$, where ${\mathcal{F}}_k={\mathcal{F}}_k^{\left(1\right)}\times {\mathcal{F}}_k^{\left(2\right)}\subset \mathcal{F}$, $k=0, 1, 2,\dots ,T$ is a corresponding increasing filtration;

• $P\left({𝓌}^{\left(1\right)},{𝓌}^{\left(2\right)}\right)=P^{\left(1\right)}\left({𝓌}^{\left(1\right)}\right)P^{\left(2\right)}\left({𝓌}^{\left(2\right)}\right)$ represents the probability measure on $\mathcal{F}$, where ${𝓌}^{\left(1\right)}\in {\mathsf{\Omega }}^{\left(1\right)}$ and ${𝓌}^{\left(2\right)}\in {\mathsf{\Omega }}^{\left(2\right)}$ are generic states of the world representing the states of the financial market and earthquake risk variables.

To explain the ZCEB price model in Equation (11), a brief explanation of the steps developed by Cox and Pedersen [55] is shown first. In the context of the comprehensive model that relies on financial market variables or earthquake risk variables, the new increasing filtrations ${\mathcal{A}}_t^{\left(1\right)}\subset {\mathcal{A}}^{\left(1\right)}$ and ${\mathcal{A}}_t^{\left(2\right)}\subset {\mathcal{A}}^{\left(2\right)}$ are defined for $t\in \left[0,T\right]$. These filtrations are generated from sigma field ${\mathcal{A}}_t^{\left(1\right)}={\mathcal{F}}_t^{\left(1\right)}\times \{\varnothing ,\mathsf{\Omega }\}$ and ${\mathcal{A}}_t^{\left(2\right)}=\left\{\varnothing ,\mathsf{\Omega }\right\}\times {\mathcal{F}}_t^{\left(2\right)}$. They prove that sigma fields ${\mathcal{A}}_T^{\left(1\right)}$ and ${\mathcal{A}}_T^{\left(2\right)}$ are independent under the probability measure $P$ so that $P\left({\alpha }_1\cap {\alpha }_2\right)=P\left({\alpha }_1\right)P\left({\alpha }_2\right)$ with ${\alpha }_1=A_1\times {\mathsf{\Omega }}^{\left(2\right)}\in {\mathcal{A}}_T^{\left(1\right)}$ for some $A_1\in {\mathcal{F}}_T^{\left(1\right)}$ and ${\alpha }_2={\mathsf{\Omega }}^{\left(1\right)}\times A_2\in {\mathcal{A}}_T^{\left(2\right)}$ for some $A_2\in {\mathcal{F}}_T^{\left(2\right)}$.

A random variable $X$ on $(\mathsf{\Omega },\mathcal{F},P)$ is dependent only on financial market risk variables if $X$ is measurable with respect to ${\mathcal{A}}_T^{\left(1\right)}$. In other words, $X\left({𝓌}^{\left(1\right)},{𝓌}^{\left(2\right)}\right)=X\left({𝓌}^{\left(1\right)}\right)$. Similarly, a random variable $X$ on $(\mathsf{\Omega },\mathcal{F},P)$ is dependent only on earthquake risk variables if $X$ is measurable with respect to ${\mathcal{A}}_T^{\left(2\right)}$. In other words, $X\left({𝓌}^{\left(1\right)},{𝓌}^{\left(2\right)}\right)=X\left({𝓌}^{\left(2\right)}\right)$.

A stochastic process $M$ is evolved through dependence only on financial market risk variables if adapted to ${\mathcal{A}}^{\left(1\right)}$. Similarly, a stochastic process $M$ is evolved through dependence only on earthquake risk variables if adapted to ${\mathcal{A}}^{\left(2\right)}$.

They proved that under the assumption that aggregate consumption depends only on the financial variables, for any random variable $X$ that depends only on catastrophic risk variables, $E^Q\left(X\right)=E^P\left(X\right)$. Put simply, the distributions of the process representing the maximum strength of extreme earthquakes remain unchanged when transitioning from the risk-neutral measure $Q$ to the physical probability measure $P$. The sigma fields ${\mathcal{A}}_T^{\left(1\right)}$ and ${\mathcal{A}}_T^{\left(2\right)}$ exhibit independence under the probability measure $Q$.

Based on the primary findings of Cox and Pedersen [55], we can infer that in the context of risk-neutral pricing measure $Q$, events that only depend on financial market risk variables are independent of those that depend on earthquake risk variables. Thus, Equation (11) can be obtained as follows: $$\begin{array}{l}{V}_t=E^Q\left(e^{-{\int }_t^T{r}_{s}ds}{R}_T|{\mathcal{F}}_t\right)\\ =E^Q\left(e^{{\int }_t^T{r}_{s}ds}|{\mathcal{F}}_t\right)E^Q\left(R_T|{\mathcal{F}}_t\right),\\ =B_{CIR}\left(t,T\right)E^P\left(R_T|{\mathcal{F}}_t\right),\\ =B_{CIR}(t,T)E^P({\eta }_{1}K{1}_{0\le M_{\left[t,T\right]}<{𝒹}_1}+{\eta }_{2}K{1}_{{𝒹}_1\le M_{\left[t,T\right]}<{𝒹}_2}+\cdots +{\eta }_{\mathcal{S}}K{1}_{M_{\left[t,T\right]}\ge {𝒹}_{\mathcal{S}-1}}|{\mathcal{F}}_t),\\ =B_{CIR}\left(t,T\right)[{\eta }_{1}KP\left(0\le M_{\left[t,T\right]}<{𝒹}_1\right)+{\eta }_{2}KP\left({𝒹}_1\le M_{\left[t,T\right]}<{𝒹}_2\right)+\cdots \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\eta }_{\mathcal{S}}KP\left(M_{\left[t,T\right]}\ge {𝒹}_{\mathcal{S}-1}\right)],\\ =B_{CIR}\left(t,T\right)K\left\{{\eta }_1\left[P\left(M_{\left[t,T\right]}\le {𝒹}_1\right)-P\left(M_{\left[t,T\right]}\le 0\right)\right]\right.\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\eta }_2\left[P\left(M_{\left[t,T\right]}\le {𝒹}_2\right)-P\left(M_{\left[t,T\right]}\le {𝒹}_1\right)\right]+\cdots \\ \left.\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\eta }_{\mathcal{S}}\left[1-P\left(M_{\left[t,T\right]}\le {𝒹}_{\mathcal{S}-1}\right)\right]\right\},\\ =B_{CIR}\left(t,T\right)K\left[{\eta }_{1}P\left(M_{\left[t,T\right]}\le {𝒹}_1\right)-{\eta }_{1}P\left(M_{\left[t,T\right]}\le 0\right)+{\eta }_{2}P\left(M_{\left[t,T\right]}\le {𝒹}_2\right)\right.\\ \left.\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-{\eta }_{2}P\left(M_{\left[t,T\right]}\le {𝒹}_1\right)+\cdots +{\eta }_{\mathcal{S}}-{\eta }_{\mathcal{S}}P\left(M_{\left[t,T\right]}\le {𝒹}_{\mathcal{S}-1}\right)\right],\\ =B_{CIR}\left(t,T\right)K[-{\eta }_{1}P\left(M_{\left[t,T\right]}\le 0\right)+{\eta }_{1}P\left(M_{\left[t,T\right]}\le {𝒹}_1\right)-{\eta }_{2}P\left(M_{\left[t,T\right]}\le {𝒹}_1\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\eta }_{2}P\left(M_{\left[t,T\right]}\le {𝒹}_2\right)-{\eta }_{3}P\left(M_{\left[t,T\right]}\le {𝒹}_2\right)+\cdots \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\eta }_{\mathcal{S}-1}P\left(M_{\left[t,T\right]}\le {𝒹}_{\mathcal{S}-1}\right)-{\eta }_{\mathcal{S}}P\left(M_{\left[t,T\right]}\le {𝒹}_{\mathcal{S}-1}\right)+{\eta }_{\mathcal{S}}],\\ =B_{CIR}\left(t,T\right)K[-{\eta }_{1}P\left(M_{\left[t,T\right]}\le 0\right)+({\eta }_1-{\eta }_2)P\left(M_{\left[t,T\right]}\le {𝒹}_1\right)+({\eta }_2\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-{\eta }_3)P\left(M_{\left[t,T\right]}\le {𝒹}_2\right)+\cdots +\left({\eta }_{\mathcal{S}-1}-{\eta }_{\mathcal{S}}\right)P\left(M_{\left[t,T\right]}\le {𝒹}_{\mathcal{S}-1}\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\eta }_{\mathcal{S}}],\\ =B_{CIR}\left(t,T\right)K[({\eta }_1-{\eta }_2)P\left(M_{\left[t,T\right]}\le {𝒹}_1\right)+({\eta }_2-{\eta }_3)P\left(M_{\left[t,T\right]}\le {𝒹}_2\right)+\cdots \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\left({\eta }_{\mathcal{S}-1}-{\eta }_{\mathcal{S}}\right)P\left(M_{\left[t,T\right]}\le {𝒹}_{\mathcal{S}-1}\right)+{\eta }_{\mathcal{S}}],\\ =B_{CIR}\left(t,T\right)K\left[{\eta }_{\mathcal{S}}+{\displaystyle \sum _{𝓂=1}^{\mathcal{S}-1}}({\eta }_{𝓂}-{\eta }_{𝓂+1})P({𝓂}_{\left[t,T\right]}\le {𝒹}_{𝓂})\right],\end{array}$$ where $$P\left(M_{\left[t,T\right]}\le {𝒹}_{𝓂}\right)=\left\{\begin{array}{cc}{e}^{-m_{\left[t,T\right]}{\left(1+\frac{\xi }{{\sigma }_2}\left({𝒹}_{𝓂}-\mu \right)\right)}^{-\frac{1}{\xi }}}& :\xi \ne 0\\ e^{-m_{\left[t,T\right]}{e}^{\left(-\frac{{𝒹}_{𝓂}-\mu }{{\sigma }_2}\right)}}& :\xi =0\end{array}\right..$$ and $$\begin{array}{l}{B}_{CIR}\left(t,T\right)=A\left(t,T\right)e^{-B\left(t,T\right)r_t},\\ A\left(t,T\right)={\left[\frac{2\gamma e^{\frac{\left(\kappa +\gamma \right)\left(T-t\right)}{2}}}{\left(\kappa +\gamma \right)\left(e^{\gamma \left(T-t\right)}-1\right)+2\gamma }\right]}^{\frac{2\kappa \theta }{{\sigma }_1^2}},\\ B\left(t,T\right)=\frac{2\left(e^{\gamma \left(T-t\right)}-1\right)}{\left(\kappa +\gamma \right)\left(e^{\gamma \left(T-t\right)}-1\right)+2\gamma },\\ \gamma =\sqrt{{\kappa }^2+2{\sigma }_1^2}.\end{array}$$

Appendix C

$$\begin{array}{l}{V}_t=E^Q\left(e^{-{\int }_t^T{r}_{s}ds}{R}_T|{\mathcal{F}}_t\right)\\ =E^Q\left(e^{{\int }_t^T{r}_{s}ds}|{\mathcal{F}}_t\right)E^Q\left(R_T|{\mathcal{F}}_t\right),\\ =B_{CIR}\left(t,T\right)E^P\left(R_T|{\mathcal{F}}_t\right),\\ =B_{CIR}(t,T)E^P({\eta }_1(K+C)1_{0\le M_{\left[t,T\right]}<{𝒹}_1}+{\eta }_2(K+C)1_{{𝒹}_1\le M_{\left[t,T\right]}<{𝒹}_2}+\cdots +{\eta }_{\mathcal{S}}(K\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+C)1_{M_{\left[t,T\right]}\ge {𝒹}_{\mathcal{S}-1}}|{\mathcal{F}}_t),\\ =B_{CIR}\left(t,T\right)[{\eta }_1(K+C)P\left(0\le M_{\left[t,T\right]}<{𝒹}_1\right)+{\eta }_2(K+C)P\left({𝒹}_1\le M_{\left[t,T\right]}<{𝒹}_2\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\cdots +{\eta }_{\mathcal{S}}(K+C)P\left(M_{\left[t,T\right]}\ge {𝒹}_{\mathcal{S}-1}\right)],\\ =B_{CIR}\left(t,T\right)(K+C)\{{\eta }_1\left[P\left(M_{\left[t,T\right]}\le {𝒹}_1\right)-P\left(M_{\left[t,T\right]}\le 0\right)\right]\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\eta }_2\left[P\left(M_{\left[t,T\right]}\le {𝒹}_2\right)-P\left(M_{\left[t,T\right]}\le {𝒹}_1\right)\right]+\cdots \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\eta }_{\mathcal{S}}\left[1-P\left(M_{\left[t,T\right]}\le {𝒹}_{\mathcal{S}-1}\right)\right]\},\\ =B_{CIR}\left(t,T\right)(K+C)[{\eta }_{1}P\left(M_{\left[t,T\right]}\le {𝒹}_1\right)-{\eta }_{1}P\left(M_{\left[t,T\right]}\le 0\right)+{\eta }_{2}P\left(M_{\left[t,T\right]}\le {𝒹}_2\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-{\eta }_{2}P\left(M_{\left[t,T\right]}\le {𝒹}_1\right)+\cdots +{\eta }_{\mathcal{S}}-{\eta }_{\mathcal{S}}P\left(M_{\left[t,T\right]}\le {𝒹}_{\mathcal{S}-1}\right)],\\ =B_{CIR}\left(t,T\right)(K+C)[-{\eta }_{1}P\left(M_{\left[t,T\right]}\le 0\right)+{\eta }_{1}P\left(M_{\left[t,T\right]}\le {𝒹}_1\right)-{\eta }_{2}P\left(M_{\left[t,T\right]}\le {𝒹}_1\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\eta }_{2}P\left(M_{\left[t,T\right]}\le {𝒹}_2\right)-{\eta }_{3}P\left(M_{\left[t,T\right]}\le {𝒹}_2\right)+\cdots \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\eta }_{\mathcal{S}-1}P\left(M_{\left[t,T\right]}\le {𝒹}_{\mathcal{S}-1}\right)-{\eta }_{\mathcal{S}}P\left(M_{\left[t,T\right]}\le {𝒹}_{\mathcal{S}-1}\right)+{\eta }_{\mathcal{S}}],\\ =B_{CIR}\left(t,T\right)(K+C)[-{\eta }_{1}P\left(M_{\left[t,T\right]}\le 0\right)+({\eta }_1-{\eta }_2)P\left(M_{\left[t,T\right]}\le {𝒹}_1\right)+({\eta }_2\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-{\eta }_3)P\left(M_{\left[t,T\right]}\le {𝒹}_2\right)+\cdots +\left({\eta }_{\mathcal{S}-1}-{\eta }_{\mathcal{S}}\right)P\left(M_{\left[t,T\right]}\le {𝒹}_{\mathcal{S}-1}\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\eta }_{\mathcal{S}}],\\ =B_{CIR}\left(t,T\right)(K+C)[({\eta }_1-{\eta }_2)P\left(M_{\left[t,T\right]}\le {𝒹}_1\right)+({\eta }_2-{\eta }_3)P\left(M_{\left[t,T\right]}\le {𝒹}_2\right)+\cdots \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\left({\eta }_{\mathcal{S}-1}-{\eta }_{\mathcal{S}}\right)P\left({𝓂}_{\left[t,T\right]}\le {𝒹}_{\mathcal{S}-1}\right)+{\eta }_{\mathcal{S}}],\\ =B_{CIR}\left(t,T\right)(K+C)\left[{\eta }_{\mathcal{S}}+{\displaystyle \sum _{𝓂=1}^{\mathcal{S}-1}}({\eta }_{𝓂}-{\eta }_{𝓂+1})P({𝓂}_{\left[t,T\right]}\le {𝒹}_{𝓂})\right],\end{array}$$ where $$P\left(M_{\left[t,T\right]}\le {𝒹}_{𝓂}\right)=\left\{\begin{array}{cc}{e}^{-m_{\left[t,T\right]}{\left(1+\frac{\xi }{{\sigma }_2}\left({𝒹}_{𝓂}-\mu \right)\right)}^{-\frac{1}{\xi }}}& :\xi \ne 0\\ e^{-m_{\left[t,T\right]}{e}^{\left(-\frac{{𝒹}_{𝓂}-\mu }{{\sigma }_2}\right)}}& :\xi =0\end{array}\right..$$ and $$\begin{array}{l}{B}_{CIR}\left(t,T\right)=A\left(t,T\right)e^{-B\left(t,T\right)r_t},\\ A\left(t,T\right)={\left[\frac{2\gamma e^{\frac{\left(\kappa +\gamma \right)\left(T-t\right)}{2}}}{\left(\kappa +\gamma \right)\left(e^{\gamma \left(T-t\right)}-1\right)+2\gamma }\right]}^{\frac{2\kappa \theta }{{\sigma }_1^2}},\\ B\left(t,T\right)=\frac{2\left(e^{\gamma \left(T-t\right)}-1\right)}{\left(\kappa +\gamma \right)\left(e^{\gamma \left(T-t\right)}-1\right)+2\gamma },\\ \gamma =\sqrt{{\kappa }^2+2{\sigma }_1^2}.\end{array}$$

References

1. Rashid, M.; Ahmad, N. Economic Losses Due to Earthquake—Induced Structural Damages in RC SMRF Structures. Cogent Eng.; 2017; 4, 1296529. [DOI: https://dx.doi.org/10.1080/23311916.2017.1296529]

2. Shin, J.Y.; Chen, S.; Kim, T.-W. Application of Bayesian Markov Chain Monte Carlo Method with Mixed Gumbel Distribution to Estimate Extreme Magnitude of Tsunamigenic Earthquake. KSCE J. Civ. Eng.; 2015; 19, pp. 366-375. [DOI: https://dx.doi.org/10.1007/s12205-015-0430-0]

3. Podolskiy, E.A.; Nishimura, K.; Abe, O.; Chernous, P.A. Earthquake-Induced Snow Avalanches: I. Historical Case Studies. J. Glaciol.; 2010; 56, pp. 431-446. [DOI: https://dx.doi.org/10.3189/002214310792447815]

4. Kalfin,; Sukono,; Supian, S.; Mamat, M. Insurance as an Alternative for Sustainable Economic Recovery after Natural Disasters: A Systematic Literature Review. Sustainability; 2022; 14, 4349. [DOI: https://dx.doi.org/10.3390/su14074349]

5. Lee, J.-P.; Yu, M.-T. Valuation of Catastrophe Reinsurance with Catastrophe Bonds. Insur. Math. Econ.; 2007; 41, pp. 264-278. [DOI: https://dx.doi.org/10.1016/j.insmatheco.2006.11.003]

6. Ismail, E. The Complementary Compound Truncated Poisson-Weibull Distribution for Pricing Catastrophic Bonds for Extreme Earthquakes. Br. J. Econ. Manag. Trade; 2016; 14, pp. 1-9. [DOI: https://dx.doi.org/10.9734/BJEMT/2016/26775] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27240745]

7. Wu, D.; Zhou, Y. Catastrophe Bond and Risk Modeling: A Review and Calibration Using Chinese Earthquake Loss Data. Hum. Ecol. Risk Assess. Int. J.; 2010; 16, pp. 510-523. [DOI: https://dx.doi.org/10.1080/10807031003788592]

8. Zimbidis, A.A.; Frangos, N.E.; Pantelous, A.A. Modeling Earthquake Risk via Extreme Value Theory and Pricing the Respective Catastrophe Bonds. ASTIN Bull.; 2007; 37, pp. 163-183. [DOI: https://dx.doi.org/10.2143/AST.37.1.2020804]

9. Kiohos, A.; Paspati, M. Alternative to Insurance Risk Transfer: Creating a Catastrophe Bond for Romanian Earthquakes. Bull. Appl. Econ.; 2021; 8, pp. 1-17. [DOI: https://dx.doi.org/10.47260/bae/811]

10. Wei, L.; Liu, L.; Hou, J. Pricing Hybrid-Triggered Catastrophe Bonds Based on Copula-EVT Model. Quant. Financ. Econ.; 2022; 6, pp. 223-243. [DOI: https://dx.doi.org/10.3934/QFE.2022010]

11. Grosser, H.; Baumbach, M.; Berckhemer, H.; Baier, B.; Karahan, A.; Schelle, H.; Krüger, F.; Paulat, A.; Michel, G.; Demirtas, R. et al. The Erzincan (Turkey) Earthquake (Ms 6.8) of March 13, 1992 and Its Aftershock Sequence. Pure Appl. Geophys.; 1998; 152, pp. 465-505. [DOI: https://dx.doi.org/10.1007/s000240050163]

12. Ioki, K.; Tanioka, Y.; Kawakami, G.; Kase, Y.; Nishina, K.; Hirose, W.; Hayashi, K.; Takahashi, R. Fault Model of the 12th Century Southwestern Hokkaido Earthquake Estimated from Tsunami Deposit Distributions. Earth Planets Space; 2019; 71, 54. [DOI: https://dx.doi.org/10.1186/s40623-019-1034-6]

13. Tang, Q.; Yuan, Z. CAT Bond Pricing under a Product Probability Measure with POT Risk Characterization. ASTIN Bull.; 2019; 49, pp. 457-490. [DOI: https://dx.doi.org/10.1017/asb.2019.11]

14. Romaniuk, M. Pricing the Risk-Transfer Financial Instruments via Monte Carlo Methods. Syst. Anal. Model. Simul.; 2003; 43, pp. 1043-1064. [DOI: https://dx.doi.org/10.1080/0232929031000150382]

15. Ibrahim, R.A.; Sukono,; Napitupulu, H.; Ibrahim, R.I. How to Price Catastrophe Bonds for Sustainable Earthquake Funding? A Systematic Review of the Pricing Framework. Sustainability; 2023; 15, 7705. [DOI: https://dx.doi.org/10.3390/su15097705]

16. Ibrahim, R.A.; Sukono,; Napitupulu, H.; Ibrahim, R.I.; Johansyah, M.D.; Saputra, J. Estimating Flood Catastrophe Bond Prices Using Approximation Method of the Loss Aggregate Distribution: Evidence from Indonesia. Decis. Sci. Lett.; 2023; 12, pp. 179-190. [DOI: https://dx.doi.org/10.5267/j.dsl.2023.3.001]

17. Sukono,; Ibrahim, R.A.; Saputra, M.P.A.; Hidayat, Y.; Juahir, H.; Prihanto, I.G.; Halim, N.B.A. Modeling Multiple-Event Catastrophe Bond Prices Involving the Trigger Event Correlation, Interest, and Inflation Rates. Mathematics; 2022; 10, 4685. [DOI: https://dx.doi.org/10.3390/math10244685]

18. Härdle, W.K.; Cabrera, B.L. Calibrating CAT Bonds for Mexican Earthquakes. J. Risk Insur.; 2010; 77, pp. 625-650. [DOI: https://dx.doi.org/10.1111/j.1539-6975.2010.01355.x]

19. Grigoryan, A. Sovereign Parametric Catastrophe Bonds as an Example of Disaster Risk Financing Instrument in the Region of ECIS State of Play and Opportunities; United Nations Development Programme: Istanbul, Türkiye, 2020.

20. Kunreuther, H.; Heal, G. Managing Catastrophic Risk. Encyclopedia of Energy, Natural Resources and Environmental Economics; Elsevier: Amsterdam, The Netherlands, 2013; pp. 52-59.

21. Deng, G.; Liu, S.; Li, L.; Deng, C. Research on the Pricing of Global Drought Catastrophe Bonds. Math. Probl. Eng.; 2020; 2020, 3898191. [DOI: https://dx.doi.org/10.1155/2020/3898191]

22. Ando, S.; Fu, C.; Roch, F.; Wiriadinata, U. Sovereign Climate Debt Instruments: An Overview of the Green and Catastrophe Bond Markets. Staff Climate Notes; IMF eLibrary: Washington, DC, USA, 2022; Volume 2022, pp. 1-28. [DOI: https://dx.doi.org/10.5089/9798400210006.066]

23. Goda, K. Multi-Hazard Parametric Catastrophe Bond Trigger Design for Subduction Earthquakes and Tsunamis. Earthq. Spectra; 2021; 37, pp. 1827-1848. [DOI: https://dx.doi.org/10.1177/8755293020981974]

24. Anggraeni, W.; Supian, S.; Sukono,; Halim, N.B.A. Earthquake Catastrophe Bond Pricing Using Extreme Value Theory: A Mini-Review Approach. Mathematics; 2022; 10, 4196. [DOI: https://dx.doi.org/10.3390/math10224196]

25. Tao, Z.; Tao, X.; Li, P. Pricing Model for Earthquake CAT Bonds. Proceedings of the 2009 International Conference on Business Intelligence and Financial Engineering; Beijing, China, 24–26 July 2009; IEEE: Piscataway, NJ, USA, 2009; pp. 740-744.

26. Jarrow, R.A. A Simple Robust Model for Cat Bond Valuation. Financ. Res. Lett.; 2010; 7, pp. 72-79. [DOI: https://dx.doi.org/10.1016/j.frl.2010.02.005]

27. Nowak, P.; Romaniuk, M. Pricing and Simulations of Catastrophe Bonds. Insur. Math. Econ.; 2013; 52, pp. 18-28. [DOI: https://dx.doi.org/10.1016/j.insmatheco.2012.10.006]

28. Ma, Z.-G.; Ma, C.-Q. Pricing Catastrophe Risk Bonds: A Mixed Approximation Method. Insur. Math. Econ.; 2013; 52, pp. 243-254. [DOI: https://dx.doi.org/10.1016/j.insmatheco.2012.12.007]

29. Chaubey, Y.P.; Garrido, J.; Trudeau, S. On the Computation of Aggregate Claims Distributions: Some New Approximations. Insur. Math. Econ.; 1998; 23, pp. 215-230. [DOI: https://dx.doi.org/10.1016/S0167-6687(98)00029-8]

30. Liu, J.; Xiao, J.; Yan, L.; Wen, F. Valuing Catastrophe Bonds Involving Credit Risks. Math. Probl. Eng.; 2014; 2014, 563086. [DOI: https://dx.doi.org/10.1155/2014/563086]

31. Jarrow, R.A.; Turnbull, S.M. Pricing Derivatives on Financial Securities Subject to Credit Risk. J. Financ.; 1995; 50, pp. 53-85. [DOI: https://dx.doi.org/10.1111/j.1540-6261.1995.tb05167.x]

32. Ma, Z.; Ma, C.; Xiao, S. Pricing Zero-Coupon Catastrophe Bonds Using EVT with Doubly Stochastic Poisson Arrivals. Discret. Dyn. Nat. Soc.; 2017; 2017, 3279647. [DOI: https://dx.doi.org/10.1155/2017/3279647]

33. Black, F.; Derman, E.; Toy, W. A One-Factor Model of Interest Rates and Its Application to Treasury Bond Options. Financ. Anal. J.; 1990; 46, pp. 33-39. [DOI: https://dx.doi.org/10.2469/faj.v46.n1.33]

34. Nowak, P.; Romaniuk, M. Valuing Catastrophe Bonds Involving Correlation and CIR Interest Rate Model. Comput. Appl. Math.; 2018; 37, pp. 365-394. [DOI: https://dx.doi.org/10.1007/s40314-016-0348-2]

35. Chao, W.; Zou, H. Multiple-Event Catastrophe Bond Pricing Based on CIR-Copula-POT Model. Discret. Dyn. Nat. Soc.; 2018; 2018, 5068480. [DOI: https://dx.doi.org/10.1155/2018/5068480]

36. Shao, J.; Pantelous, A.; Papaioannou, A.D. Catastrophe Risk Bonds with Applications to Earthquakes. Eur. Actuar. J.; 2015; 5, pp. 113-138. [DOI: https://dx.doi.org/10.1007/s13385-015-0104-9]

37. Gunardi,; Setiawan, E.P. Valuation of Indonesian Catastrophic Earthquake Bonds with Generalized Extreme Value (GEV) Distribution and Cox-Ingersoll-Ross (CIR) Interest Rate Model. Proceedings of the 2014 International Conference on Actuarial Science and Statistics, AIP Conference Proceedings; Bandung, Indonesia, 21–23 October 2014; AIP Publishing: Melville, NY, USA, 2015; pp. 1-14.

38. Hofer, L.; Zanini, M.A.; Gardoni, P. Risk-Based Catastrophe Bond Design for a Spatially Distributed Portfolio. Struct. Saf.; 2020; 83, 101908. [DOI: https://dx.doi.org/10.1016/j.strusafe.2019.101908]

39. Mistry, H.K.; Lombardi, D. Pricing Risk-Based Catastrophe Bonds for Earthquakes at an Urban Scale. Sci. Rep.; 2022; 12, 9729. [DOI: https://dx.doi.org/10.1038/s41598-022-13588-1]

40. Kang, H.-B.; Chang, H.; Chang, T. Catastrophe Reinsurance Pricing-Modification of Dynamic Asset-Liability Management. Rom. J. Econ. Forecast; 2022; 25, pp. 5-20.

41. Aghdam, Y.E.; Neisy, A.; Adl, A. Simulating and Pricing CAT Bonds Using the Spectral Method Based on Chebyshev Basis. Comput. Econ.; 2022; 63, pp. 423-435. [DOI: https://dx.doi.org/10.1007/s10614-022-10347-2]

42. Anggraeni, W.; Supian, S.; Sukono,; Halim, N.A. Single Earthquake Bond Pricing Framework with Double Trigger Parameters Based on Multi Regional Seismic Information. Mathematics; 2023; 11, 689. [DOI: https://dx.doi.org/10.3390/math11030689]

43. Hofer, L.; Gardoni, P.; Zanini, M.A. Risk-Based CAT Bond Pricing Considering Parameter Uncertainties. Sustain. Resilient Infrastruct.; 2021; 6, pp. 315-329. [DOI: https://dx.doi.org/10.1080/23789689.2019.1667116]

44. Cummins, J.D. CAT Bonds and Other Risk-Linked Securities: State of the Market and Recent Developments. SSRN Electron. J.; 2007; 11, pp. 23-47. [DOI: https://dx.doi.org/10.2139/ssrn.1057401]

45. Loubergé, H.; Kellezi, E.; Gilli, M. Using Catastrophe-Linked Securities to Diversify Insurance Risk: A Financial Analysis of Cat Bonds. J. Insur. Issues; 1999; 22, pp. 125-146.

46. Cummins, J.D.; Weiss, M.A. Convergence of Insurance and Financial Markets: Hybrid and Securitized Risk-Transfer Solutions. J. Risk Insur.; 2009; 76, pp. 493-545. [DOI: https://dx.doi.org/10.1111/j.1539-6975.2009.01311.x]

47. Burnecki, K.; Kukla, G.; Taylor, D. Pricing of Catastrophe Bonds. Statistical Tools for Finance and Insurance; Cizek, P.; Härdle, W.; Weron, R. Springer: Berlin/Heidelberg, Germany, 2005; pp. 93-114.

48. Cox, J.C.; Ingersoll, J.E.; Ross, S.A. A Theory of the Term Structure of Interest Rates. Econometrica; 1985; 53, pp. 385-407. [DOI: https://dx.doi.org/10.2307/1911242]

49. Vasicek, O. An Equilibrium Characterization of the Term Structure. J. Financ. Econ.; 1977; 5, pp. 177-188. [DOI: https://dx.doi.org/10.1016/0304-405X(77)90016-2]

50. Salvadori, G.; De Michele, C.; Kottegoda, N.T.; Rosso, R. Univariate Extreme Value Theory. Extremes in Nature; Springer: Dordrecht, The Netherlands, 2007; Volume 56, pp. 1-112.

51. Balkema, A.A.; de Haan, L. Residual Life Time at Great Age. The Annals of Probability; Institute of Mathematical Statistics: Waite Hill, OH, USA, 1974; Volume 2, [DOI: https://dx.doi.org/10.1214/aop/1176996548]

52. Ibrahim, R.A.; Sukono,; Riaman,. Estimation of the Extreme Distribution Model of Economic Losses Due to Outbreaks Using the POT Method with Newton Raphson Iteration. Int. J. Quant. Res. Model.; 2021; 2, pp. 37-45. [DOI: https://dx.doi.org/10.46336/ijqrm.v2i1.118]

53. Osaki, S. Poisson Processes. Applied Stochastic System Modeling; Springer: Berlin/Heidelberg, Germany, 1992; pp. 63-82.

54. Klugman, S.A.; Panjer, H.H.; Willmot, G.E. Aggregate Loss Models. Loss Models: From Data to Decisions; John Wiley & Sons, Inc.: New York, NY, USA, 2008; pp. 199-268.

55. Cox, S.H.; Pedersen, H.W. Catastrophe Risk Bonds. N. Am. Actuar. J.; 2000; 4, pp. 56-82. [DOI: https://dx.doi.org/10.1080/10920277.2000.10595938]

56. Setiawan, E.P.; Wutsqa, D.U.; Abadi, A.M.; Kusuma, E. Pricing Indonesian Earthquake Catastrophe Bond Based on Depth and Magnitude; AIP Publishing: Melville, NY, USA, 2023; 080025.

57. Sukono,; Napitupulu, H.; Riaman,; Ibrahim, R.A.; Johansyah, M.D.; Hidayana, R.A. A Regional Catastrophe Bond Pricing Model and Its Application in Indonesia’s Provinces. Mathematics; 2023; 11, 3825. [DOI: https://dx.doi.org/10.3390/math11183825]

58. Zulfikar, Z.; Syahnur, S.; Majid, M.S.A. The Effect of Energy Consumption, Energy Resources, Economic Growth, and Road Infrastructure on CO₂ Emissions in Indonesia. Int. J. Quant. Res. Model.; 2021; 2, pp. 173-183. [DOI: https://dx.doi.org/10.46336/ijqrm.v2i3.173]

59. Thode, H.C. Testing For Normality; 1st ed. CRC Press: Boca Raton, FL, USA, 2002; ISBN 9780203910894

60. Anderson, T.W.; Darling, D.A. Asymptotic Theory of Certain “Goodness of Fit” Criteria Based on Stochastic Processes. Ann. Math. Stat.; 1952; 23, pp. 193-212. [DOI: https://dx.doi.org/10.1214/aoms/1177729437]

61. Dodge, Y. Chi-Square Goodness of Fit Test. The Concise Encyclopedia of Statistics; Springer: New York, NY, USA, 2008; pp. 72-76.

62. Massey, F.J. The Kolmogorov-Smirnov Test for Goodness of Fit. J. Am. Stat. Assoc.; 1951; 46, 68. [DOI: https://dx.doi.org/10.1080/01621459.1951.10500769]

63. Vere-Jones, D. On the Estimation of Frequency in Point-Process Data. J. Appl. Probab.; 1982; 19, pp. 383-394. [DOI: https://dx.doi.org/10.2307/3213577]

64. Merdasse, M.; Hamdache, M.; Peláez, J.A.; Henares, J.; Medkour, T. Earthquake Magnitude and Frequency Forecasting in Northeastern Algeria Using Time Series Analysis. Appl. Sci.; 2023; 13, 1566. [DOI: https://dx.doi.org/10.3390/app13031566]

65. Lewis, C.D. Industrial and Business Forecasting Methods; Butterworth-Heinemann: London, UK, 1982.

66. Alfonsi, A. On the Discretization Schemes for the CIR (and Bessel Squared) Processes. Monte Carlo Methods Appl.; 2005; 11, pp. 355-384. [DOI: https://dx.doi.org/10.1515/156939605777438569]