2.1. Fractional Brownian motion

We consider a filtered probability space $(\Omega, \mathcal{F}, \mathcal{F}_t, \mathbb{P})$ such that all random variables are well defined and stochastic processes are $\mathcal{F}_t$ adapted. The natural filtration of a stochastic process is taken as the $\mathbb{P}$ -completion of the filtration that it generates. For some fixed $T\gt0$ , a (standard) fractional Brownian motion (fBm) on time interval [0, T], $B^H = (B_t^H, t \in [0, T])$ with Hurst parameter $H \in (0, 1)$ is defined as a centered Gaussian process with continuous paths such that $B_0^H = 0$ and

\begin{align*}\mathrm{Cov}(B^H_s, B^H_t) = \mathbb{E}\left(B^H_s B^H_t\right) = \frac{1}{2}\left(s^{2H} + t^{2H} - |s-t|^{2H}\right), \ 0 \leq s, t \leq T. \end{align*}

The Hurst parameter $H \in (0,1)$ governs the memory structure of the fBm and has a direct economic interpretation. When $H = \tfrac{1}{2}$ , the fBm reduces to standard Brownian motion, where increments are independent and the process has no memory of past shocks. When $H \gt \tfrac{1}{2}$ , the process exhibits LRD. In the context of excess mortality, $H \gt \tfrac{1}{2}$ captures the phenomenon whereby a mortality shock persists over many periods rather than reverting quickly, such as the lasting effect of the COVID-19 pandemic. For interest rates, $H \gt \tfrac{1}{2}$ reflects prolonged monetary-policy regimes, such as the low interest rate periods observed over 2008–2015 and 2020–2021. Conversely, when $H \lt \tfrac{1}{2}$ , the process is antipersistent, meaning that positive increments tend to be followed by negative corrections. This will produce rougher and more rapidly mean-reverting trajectories that may be relevant for short-term fluctuations that are quickly corrected by market forces. Lastly, we remark that, when $H \neq \frac{1}{2}$ , fBm is not a semimartingale and the classic Itô theory cannot be applied to study the corresponding stochastic integrals.

The construction of stochastic integrals can follow the traditional limiting arguments (see Hu et al., Reference Hu, Nualart and Zhou2019; Nualart, Reference Nualart2006 and the references therein). Specifically, we consider the elementary space $\mathcal{E}$ of real valued step functions and define the Hilbert space $\mathfrak{H}$ as the closure of the space $\mathcal{E}$ endowed with the inner product

\begin{align*} \langle \unicode{x1D7D9}_{[0, t]}, \unicode{x1D7D9}_{[0, s]} \rangle_{\mathfrak{H}} = \mathbb{E}(B_s^H B_t^H). \end{align*}

The mapping $\unicode{x1D7D9}_{[0, t]} \to B_t^H$ can be extended to the linear isometry between the Hilbert space and the Gaussian subspace spanned by fBm. In other words, $\int_0^T \unicode{x1D7D9}_{[0, t]} (s) dB_s^H = B^H_t$ . By the limiting argument, the stochastic integral (also called the Wiener integral) $ \int_0^T f(t) dB_t^H$ whose integrand $f \in \mathfrak{H}$ is deterministic follows a centered Gaussian distribution. When $H\gt\frac{1}{2}$ , the variance can be computed by

(1) \begin{equation} \mathbb{E}\left[\left(\int_0^T f(t) dB_t^H\right)^2\right] = \alpha_H \int_0^T \int_0^T f(s) f(t) |t-s|^{2H-2} {\mathop{}\! {\mathop{}\!\mathrm{d}}s{\mathop{}\!\mathrm{d}}t}, \end{equation}

where $\alpha_H = H(2H-1)$ . The covariance between two Wiener integrals is computed by

(2) \begin{equation} \mathrm{Cov}\left(\int_0^T f(t) dB_t^H, \int_0^T g(t) dB_t^H\right) = \alpha_H\int_0^T \int_0^T f(t) g(s) |t-s|^{2H-2} {\mathop{}\! {\mathop{}\!\mathrm{d}}s{\mathop{}\!\mathrm{d}}t}. \end{equation}

There are some integral transformations that can change fBm to martingales (see Norros et al., Reference Norros, Valkeila and Virtamo1999). The following results from Norros et al. (Reference Norros, Valkeila and Virtamo1999) will be used in this paper. We first introduce the deterministic functions for $0 \lt s \lt t \leq T$ ,

\begin{align*}k_H(t, s) = \kappa_H^{-1} s^{\frac{1}{2} - H} (t-s)^{\frac{1}{2} - H}, \ \kappa_H = 2H\Gamma\bigg(\frac{3}{2}-H\bigg)\Gamma\bigg(H+\frac{1}{2}\bigg), \end{align*}

\begin{align*}\omega_t^H = \lambda_H^{-1} t^{2-2H}, \ \lambda_H = \frac{2H\Gamma(3-2H)\Gamma(H+\frac{1}{2})}{\Gamma(\frac{3}{2}-H)}. \end{align*}

Then, we can construct a Gaussian martingale $M_t^H = \int_0^t k_H(t, s) dB_s^H$ , and it is called the fundamental martingale. The natural filtration generated by this martingale actually coincides with that of the fBm $B^H$ , and the quadratic variation $\langle M^H \rangle_t$ equals $\omega_t^H$ . The fundamental martingale is very useful in the stochastic analysis about fBm, and one of the applications is to derive the Girsanov formula for fBm (see Kleptsyna et al., Reference Kleptsyna, Breton and Roubaud2000). Namely, we consider a process $Y_t$ that is the pathwise solution to the following stochastic differential equation

\begin{align*}dY_t = A(t) dt + \sigma(t) dB_t^H, \ 0 \leq t \leq T, \end{align*}

where the functions $A(t), \sigma(t)$ satisfy some nice properties to ensure the existence of a pathwise solution for a given initial condition. The process defined by

\begin{align*}Q_H(t) = \frac{d}{d \omega_t^H} \int_0^t k_H(t, s)\sigma^{-1}(s)A(s) {\mathop{}\! {\mathop{}\!\mathrm{d}}s}, \ t \in [0, T], \end{align*}

belongs to the measurable space $L^2([0, T], d\omega^H)$ . Then, we can construct the random variable

\begin{align*}Z_H(T) = \exp\left(-\int_0^T Q_H(t) dM_t^H - \frac{1}{2} \int_0^T (Q_H(t))^2 d\omega_t^H\right) \end{align*}

such that $\mathbb{E}[Z_H(T)] = 1$ . Consequently, $\mathbb{Q} = Z(T)\mathbb{P}$ is a probability measure, and the process Y under the measure $\mathbb{Q}$ satisfies the following equation

\begin{align*}dY_t = \sigma(t) dB_t^{H, \mathbb{Q}}, \ t \in [0, T], \end{align*}

where $B_t^{H, \mathbb{Q}}$ is a fBm with Hurst parameter H under measure $\mathbb{Q}$ .

2.2. Mixed fractional Brownian motion

We will use the combination of fBm and Bm for our proposed stochastic models, as it offers the flexibility of describing different types of driving noise. A mixed fractional Brownian motion (mfBm), $B = B^{\alpha, H} = \{B_t, t \in [0, T]\}$ with parameters $\alpha$ and H, was first introduced in Cheridito (Reference Cheridito2001). It is defined as

(3) \begin{equation} B_t = \alpha W_t + B_t^H,\end{equation}

where $\{W_t, t \in [0, T]\}$ is a standard Brownian motion, $\{B_t^H, t \in [0, T]\}$ is a fBm with Hurst parameter H, and these two processes are assumed to be independent. The parameter $\alpha \in \mathbb{R}^+$ governs the relative contribution of the standard Brownian component to the mfBm.

The two processes underlying the mfBm in Equation (3) have distinct economic roles in mortality and interest rate modeling. The standard Bm captures short-lived, memoryless shocks; for example, transitory weather-related mortality fluctuations or short-term liquidity-driven movements in the money market. By contrast, the fBm captures persistent dynamics; for example, the gradual and sustained impact of pandemics and the slow build-up and unwinding of inflation expectations that shape long-term interest rate movements. By combining both components, the mfBm offers the flexibility to control the relative contribution of short-memory and long-memory effects within a single unified setup.

An mfBm inherits many properties from the fBm. It is not a semimartingale if $H \neq \frac{1}{2}$ . When $H=\frac{1}{2}$ , the mfBm is essentially a Bm. When $H \in (\frac{1}{2}, 1)$ , the increment of mfBm exhibits long-range dependence (LRD). Namely, we can take a look at the correlation of the increments. Since the mfBm has stationary increments, we can compute

\begin{align*}\rho(k) \, :\!= \, \mathbb{E}[B_1(B_{k+1} - B_k)] = \frac{1}{2}[(k+1)^{2H} - 2k^{2H} + (k-1)^{2H}], \end{align*}

for any natural number $k \geq 1$ . Clearly, when k is large, the correlation $\rho(k)$ behaves like $k^{2H-2}$ and $\sum_k \rho(k) = \infty$ when $H \gt \frac{1}{2}$ . This shows that the mfBm has LRD when $H \gt \frac{1}{2}$ .

To model the joint dynamics of interest rate and mortality rate under the influence of both temporary and persistent effect, it is natural to consider the use of mfBm. On the one hand, standard Bm is commonly employed to capture short-term stochastic fluctuations arising from the economic and ecological environment. On the other hand, there is growing empirical evidence that mortality rates and economic data exhibit long-memory behavior, as documented in Zhou et al. (Reference Zhou, Zhou and Li2022) and references therein. The LRD reflects the persistent impact of sociodemographic, medical, and environmental factors over extended periods. Similarly, studies such as Li et al. (Reference Li, Liu, Tang and Yuan2023) have shown interest rates, under stressed conditions, exhibit sensitivity to health-related shocks. They do not necessarily imply a direct causal relationship but they reveal meaningful empirical association. As noted, the COVID-19 pandemic has reinforced this interdependence and highlighted how systemic events can induce lasting correlated shifts in both financial market and human mortality. These observations suggest a modeling framework based on mfBm, in contrast to traditional Brownian models, offers a more realistic and flexible approach for capturing the joint behavior of interest rates and mortality rates. In the next section, we introduce a general version of this framework and present its mathematical formulation.

4.1. Valuation of zero-coupon bond

Zero-coupon bonds play a crucial role in the design of mortality-linked securities by serving as instruments for discounting future cash flows. To model this, we assume that the dynamics of the interest rate follows the process given in (9). At time t, the value of a risk-free zero-coupon bond maturing at time T is given byFootnote ¹

\begin{align*}P(t, T) = \mathbb{E}^{\mathbb{Q}}\left[\exp\left(-\int_t^T r_s {\mathop{}\! {\mathop{}\!\mathrm{d}}s}\right)\Big|\mathcal{F}_t\right]. \end{align*}

To calculate the conditional expectation, we need an important result that provides an explicit formula for conditional expectation for a functional of fBm. This result is recalled in the following Lemma 4.1. For any appropriate function f, we first introduce an operator $(\!-\!\Delta)^\alpha$ defined as

\begin{align*}(\!-\!\Delta)^\alpha f(x) \, :\!= \, \frac{1}{2\Gamma(-2\alpha)\cos(-\alpha \pi)} \int_{-\infty}^\infty |x-t|^{-2\alpha-1} f(t) {\mathop{}\! {\mathop{}\!\mathrm{d}}t}, \end{align*}

for $\alpha\lt0$ , and the norm of f on the Hilbert space $\mathfrak{H}$ is denoted by

\begin{align*}\|f\|_{\mathfrak{H}}^2 \, :\!= \, \alpha_H \int_0^\infty \int_0^\infty f(u)f(v)|u-v|^{2H-2} {\mathop{}\! {\mathop{}\!\mathrm{d}}u}{\mathop{}\! {\mathop{}\!\mathrm{d}}v}, \end{align*}

where $\alpha_H = H(2H-1)$ and $H \in (\frac{1}{2}, 1)$ . According to the results given in Gripenberg and Norros (Reference Gripenberg and Norros1996) and Theorem 3.2 in Hu (Reference Hu2001), the conditional expectation of a functional fBm, with respect to the filtration generated by the fBm, can be represented as a stochastic integral with respect to the same fBm.

Lemma 4.1. Suppose a given function g satisfies necessary conditions such that the stochastic integral $\int_0^\infty g(s) dB^H(s)$ is well defined, where $B^H$ is a fBm with Hurst parameter H. Then there exists a function h with support on [0, t], such that

\begin{equation*} \mathbb{E}\left[\int_0^\infty g(s) dB^H(s) \Big| \mathcal{F}_t\right] = \int_0^t h(s) dB^H(s), \end{equation*}

and the function h satisfies

(13) \begin{equation} (\!-\!\Delta)^{\frac{1}{2}-H} g(x) = (\!-\!\Delta)^{\frac{1}{2}-H} h(x), \end{equation}

for any $x \in (0, t)$ . Moreover,

\begin{align*}\mathbb{E}\left[\exp\left(\int_0^\infty g(s) dB^H(s) \right)\Big| \mathcal{F}_t\right] = \exp\left(\int_0^t h(s) dB^H(s) - \frac{1}{2}\|h_{[0, t]}\|_{\mathfrak{H}}^2 + \frac{1}{2}\|g\|^2_{\mathfrak{H}}\right). \end{align*}

Applying the above result, we can derive an explicit formula for the bond values. This is stated in the following theorem.

Theorem 4.2. Consider a risk-free zero-coupon bond issued at time 0 and matured at time T, whose interest rate $r_t$ under measure $\mathbb{Q}$ is described by the following SDE driven by a mfBm $B_t = \alpha W_t + B_t^H$ :

\begin{align*}dr_t = (m - \theta r_t) dt + \sigma dB_t, \ t \in [0, T]. \end{align*}

The bond value at time 0 can be determined by

(14) \begin{eqnarray} P(0, T) & = & \exp\left(-\frac{m}{\theta}T + \left(\frac{m}{\theta} - r_0\right)\frac{1-e^{-\theta T}}{\theta}\right) \times \exp \left(\frac{\alpha^2 \sigma^2}{2\theta^2} \left(-\frac{1}{2\theta} e^{-2\theta T} + \frac{2}{\theta}e^{-\theta T} + T - \frac{3}{2\theta}\right)\right) \notag\\[3pt] & & \times \exp \Big(\frac{\sigma^2 \alpha_H}{\theta^2} \Big(-\frac{e^{-2\theta T}}{2\theta} \int_0^T e^{\theta x} x^{2H-2} dx + \frac{1}{2\theta}\int_0^T e^{-\theta x} x^{2H-2} dx \notag\\[3pt] & & \ - \frac{1}{2H-1} \int_0^T e^{-\theta x} x^{2H-1} dx - \frac{e^{-\theta T}}{2H-1} \int_0^T e^{\theta x}x^{2H-1} dx + \frac{T^{2H}}{(2H)(2H-1)} \Big) \Big). \end{eqnarray}

The bond value at time t can be computed as

(15) \begin{eqnarray} P(t, T) & = & \exp\left(\int_0^t r_s {\mathop{}\! {\mathop{}\!\mathrm{d}}s} -\frac{mT}{\theta} + \frac{1}{\theta}\left(\frac{m}{\theta} - r_0\right)\left(1 - e^{-\theta T}\right)\right) \times \exp\left(-\frac{\alpha \sigma}{\theta} \int_0^t 1 - e^{-\theta(T-u)} dW_u\right) \notag\\[3pt] && \times \exp\left(\frac{\alpha^2 \sigma^2}{2\theta^2} \left((T-t) - \frac{e^{-2\theta(T-t)} - 4 e^{-\theta(T-t)} + 3}{2\theta}\right)\right) \notag\\[3pt] && \times \exp\left(\int_0^t h(s) dB_s^H - \frac{1}{2}\|h_{[0, t]}\|_{\mathfrak{H}}^2 + \frac{1}{2}\|g\|^2_{\mathfrak{H}}\right), \end{eqnarray}

where the function h satisfies the Equation (13) associated with the function g that is defined by $g(u) = -\frac{\sigma}{\theta}\unicode{x1D7D9}_{[0, T]} \left(1 - e^{-\theta(T-u)}\right).$

The proof of this theorem is omitted for reasons of space. Readers are referred to the full version available online. It primarily relies on the stochastic Fubini Theorem, Itô isometry and the variance formula (1).

We provide a numerical illustration in Figure 1 to show how the value of a zero-coupon bond at time 0 is influenced by different model parameters. Figure 1 shows that the bond value does not exhibit a clear exponential decay as the bond duration increases, since the interest rate is modeled by a mean-reverting process instead of a constant. Under different volatility levels ( $\sigma = 0.5$ in the left panel and $\sigma = 1$ in the right panel), the effect of the Hurst parameter H on bond values is consistent across bond durations. In general, a larger Hurst parameter leads to a higher bond value. Moreover, comparing the two panels, we can see that as the volatility parameter $\sigma$ increases, the influence of the Hurst parameter on the bond values becomes more pronounced.

Figure 1.

An illustration of zero-coupon bond value at time 0 under different values of H and $\sigma$ , where $\alpha = 0$ , $m = 10$ , and $\theta = 1$ .

4.2 Pricing mortality-linked securities

Mortality bonds are a type of mortality-linked security (MLS) that enables insurance companies and pension plans to transfer catastrophic mortality risk to capital markets. For example, Swiss Re introduced Vita Capital bonds in 2003 to hedge against increased mortality rates in the United States, United Kingdom, Canada, and Germany. The payoff structure of these bonds is tied to a specific mortality index, such as a national death rate. If this index exceeds a predetermined threshold, the principal repayment of the bond is reduced proportionally, allowing issuers to mitigate losses from unexpected mortality events such as pandemics or natural disasters. In return, investors accept this risk in exchange for yields that are typically higher than those offered by standard bonds.

Consider an MLS issued at time 0 with scheduled payments at $t_k \in \{t_1, t_2, \ldots, t_K\}$ , where $t_K = T$ is the maturity time. At each time point $t_k$ , the bond payment $\phi_{t_k}$ depends on a mortality-linked quantity. It is assumed that for all $k = 1, 2, \ldots, K$ , the payment function $\phi_{t_k}$ is measurable with respect to the augmented filtration $\mathcal{F}_{t_k}^r$ generated by the stochastic processes $\{W_{i,u},\ 0 \lt u \lt t_k\}$ and $\{B^{H_i}_{i,u},\ 0 \lt u \lt t_k\}$ for $i = 1, 2$ .

For a mortality bond structured similarly to Vita bonds, assume the face value is F and the annual coupon rate is c, payable m times per year. The payment $\phi_{t}$ at times $t \in \{t_1, t_2, \ldots, t_K\}$ is given by

(16) \begin{align}\phi_{t} = \begin{cases} \dfrac{cF}{m} & \mbox{if } t=t_1, ..., t_{K-1} \\[9pt]\dfrac{cF}{m} + F(1-\mathrm{PRF}) & \mbox{if } t=t_K \end{cases},\end{align}

where PRF is the principal reduction factor that determines the proportion of the principal repaid to investors at maturity. The PRF depends on the realized value of a mortality index $\tilde{\mu}_t$ for $t \in \{t_1, t_2, \ldots, t_K\}$ and is defined as

\begin{align*}\mathrm{PRF}(\tilde{\mu}) = \min\left(1, \; \sum_{t= t_1}^{t_K}\frac{(\tilde{\mu}_{t} - a)_+ - (\tilde{\mu}_{t} - b)_+}{b - a} \right), \end{align*}

where a is the attachment point and b is the detachment point. If $\tilde{\mu}_t \lt a$ for all t, then $\mathrm{PRF} = 0$ and the principal is fully repaid. In contrast, if $\tilde{\mu}_t \gt b$ for any t, then $\mathrm{PRF} = 1$ , meaning that the principal is entirely lost.

Regarding the mortality index linked to the bond, we present the following considerations. First, it should be designed to accurately reflect the mortality experience of the population over the term of the bond. Several specifications for $\tilde{\mu}_t$ can be considered:

• • $\tilde{\mu}_t$ is set to the mortality rate observed at time t,

• • $\tilde{\mu}_t$ is set to the average mortality rate observed between times $t-1$ and t, or

• • $\tilde{\mu}_t$ is set to the maximum mortality rate observed between times $t-1$ and t.

Second, the choice of stochastic model plays an important role in determining the mortality index, which is an essential part in the design of the MLS. In this paper, we assume that the excess mortality rate follows the SDE (4), which highlights the LRD property observed in human mortality data. Existing literature such as Zhou et al. (Reference Zhou, Zhou and Li2022) has studied how the Hurst parameter affects the distribution of lifetime related random variables. Since our focus is on mortality bonds, it is also important to understand how extreme values of mortality respond to changes in the Hurst parameter. Figure 2 illustrates the impact of Hurst parameter $H_2$ on the maximum mortality rate over the bond term. A higher value of $H_2$ leads to a higher exceedance probability when the threshold mortality rate is set high. This suggests that stronger LRD would increase the likelihood of extreme mortality rates, which further affects the design and pricing of mortality bonds.

Figure 2.

Probability that the maximum excess mortality $\tilde\mu = \max_{0 \leq t \leq T} \mu_t$ exceeds a given threshold under different values of $H_2$ . Results are based on $10^5$ simulated trajectories of $\mu_t$ from the SDE with parameters $\alpha = 0$ , $\sigma = 0.116724$ , $m = 0.008818679$ , $\theta = 1.029601$ , and $\mu_0 = -0.085646188$ .

Based on the payout structure given in Equation (16), the time-t value of future bond payments can be expressed as:

\begin{align*} P_t^{\mathrm{MB}} &= \mathbb{E}^{\mathbb{Q}}\left[ \sum_{t \lt t_k \leq T} \phi_{t_k} e^{-\int_t^{t_k} r_s\,\mathrm{d}s} \,\middle|\, \mathcal{F}_t^r \right] \\[6pt] &= \sum_{t \lt t_k \leq T} \frac{cF}{m} \mathbb{E}^{\mathbb{Q}}\left[ e^{-\int_t^{t_k} r_s\,\mathrm{d}s} \,\middle|\, \mathcal{F}_t^r \right] + \mathbb{E}^{\mathbb{Q}}\left[ F\bigl(1-\operatorname{PRF}(\tilde{\mu})\bigr) e^{-\int_t^{T} r_s\,\mathrm{d}s} \,\middle|\, \mathcal{F}_t^r \right] \\[6pt]&= \sum_{t \lt t_k \leq T}\frac{cF}{m} P(t,t_k) + F P(t,T) - F\mathbb{E}^{\mathbb{Q}}\left[ \operatorname{PRF}(\tilde{\mu}) e^{-\int_t^{T} r_s\,\mathrm{d}s} \,\middle|\, \mathcal{F}_t^r \right]. \end{align*}

where $P(t, t_k)$ , for $t \lt t_k \leq T$ , represents the time-t price of a risk-free zero-coupon bond maturing at time $t_k$ . The explicit formula for $P(t, t_k)$ has been derived in Theorem 4.2.

Finally, the price of the mortality bond at time $t = 0$ , with an annual coupon rate of c payable m times per year, is given by

(17) \begin{equation} P_{0}^{\mathrm{MB}}(c) = \frac{cF}{m} \sum_{k=1}^{K} P(0,t_k) + F P(0,T) - F \mathbb{E}^{\mathbb{Q}} \left[ \mathrm{PRF}(\tilde{\mu}) \, e^{-\int_0^T r_s\,\mathrm{d}s} \middle| \mathcal{F}_0^r \right]. \end{equation}

where the term $\mathbb{E}^{\mathbb{Q}}\left[ \mathrm{PRF}(\tilde{\mu}) \times e^{-\int_0^T r_s {\mathop{}\! {\mathop{}\!\mathrm{d}}s}} \Big| \mathcal{F}_0^r \right]$ represents the expected value of the discounted principal reduction and can be evaluated numerically using simulations.

5.1. Parameter estimation of a general SDE driven by mfBm

In this subsection, we consider the parameter estimation of a general SDE driven by mfBm,

\begin{align*}d X_t = (m - \theta X_t) dt + \sigma (\alpha dW_t + d B_t^H), \end{align*}

where $X_0 = x_0$ . Suppose that we have one observed trajectory of $X_t$ at equidistant time points $t=0, \frac{1}{n}, \frac{2}{n}, \ldots, \frac{[nT]}{n}$ . The parameters $m, \theta, \sigma, \alpha, H$ are unknown.

We employ the well-known rescaled range (R/S) analysis (see, for example, Lo, Reference Lo1991; Mandelbrot, Reference Mandelbrot1972; Mandelbrot and Wallis, Reference Mandelbrot and Wallis1969; Wang et al., Reference Wang, Xiao and Yu2023) to estimate the Hurst parameter H. Namely, for a long memory series of length N, $\left\{X_{\frac{i}{n}}, i=0, 1, \ldots, N\right\}$ , the rescaled range can be calculated as follows:

\begin{align*}Y_j = X_{\frac{j}{n}} - \frac{1}{N} \sum_{i=1}^N X_{\frac{i}{n}}, \ j=1, 2, \ldots, N; \ Z_k = \sum_{j=1}^k Y_j, k = 1, 2, \ldots, N. \end{align*}

The range R and rescaled range $R/S(N)$ can be computed by

\begin{align*}R(N) = \max(Z_1, \ldots, Z_N) - \min(Z_1, \ldots, Z_N), \ R/S(N) = \frac{R(N)}{\sqrt{\frac{1}{N}\sum_{j=1}^N Y_j^2}}. \end{align*}

Then, the Hurst parameter H can be estimated as 0.5 plus the slope of $R/S(N) - \mathbb{E}[R/S(N)]$ against N. The formula for approximating the expectation of $R/S(N)$ can be found in (3.9) of Annis and Lloyd (Reference Annis and Lloyd1976).

To estimate the volatility parameter $\sigma$ and weight parameter $\alpha$ , we define the power variations of the stochastic process $X = \{X_t, t \geq 0\}$ as

\begin{align*}V^n(X)_t = \sum_{i=1}^{[nt]}\left|\Delta X_{\frac{i-1}{n}}\right|^2 = \sum_{i=1}^{[nt]}\left|X_{\frac{i}{n}} - X_{\frac{i-1}{n}}\right|^2, \end{align*}

\begin{align*}U^n(X)_t = \sum_{i=1}^{[nt]}\left|\tilde\Delta X_{\frac{i-1}{n}}\right|^2 = \sum_{i=1}^{[nt]}\left|X_{\frac{i+1}{n}} - X_{\frac{i-1}{n}}\right|^2. \end{align*}

We may find that at time $t=T$ ,

\begin{align*}V^n(X)_T = \sigma^2 \left(\alpha^2 V^n(W)_T + V^n(B^H)_T + 2\alpha \sum_{i=1}^{[nT]}\left|\Delta W_{\frac{i-1}{n}}\right|\left|\Delta B^H_{\frac{i-1}{n}}\right|\right) + o(n^{-1}), \end{align*}

\begin{align*}U^n(X)_T = \sigma^2 \left(\alpha^2 U^n(W)_T + U^n(B^H)_T + 2\alpha \sum_{i=1}^{[nT]}\left|\tilde\Delta W_{\frac{i-1}{n}}\right|\left|\tilde\Delta B^H_{\frac{i-1}{n}}\right|\right) + o(n^{-1}). \end{align*}

From the ergodic theorem, we have the following asymptotic results (see Hu et al., Reference Hu, Nualart and Zhou2019 and the references therein),

\begin{align*} \lim_{n \to \infty} V^n(W)_T =T, \ \lim_{n \to \infty} n^{2H-1} V^n(B^H)_T = T , \end{align*}

and

\begin{align*}\lim_{n \to \infty} n^{H-\frac{1}{2}}\sum_{i=1}^{[nT]}\left|\Delta W_{\frac{i-1}{n}}\right|\left|\Delta B^H_{\frac{i-1}{n}}\right| = T, \end{align*}

in probability. Similarly,

\begin{align*} \lim_{n \to \infty} U^n(W)_T = \mathbb{E}\left[|W_2|^2\right] = 2T, \ \lim_{n \to \infty} n^{2H-1} U^n(B^H)_T = \mathbb{E}\left[\left|B_2^H\right|^2\right] = 2^{2H} T, \end{align*}

and

\begin{align*}\lim_{n \to \infty} n^{H-\frac{1}{2}}\sum_{i=1}^{[nT]}\left|\tilde\Delta W_{\frac{i-1}{n}}\right|\left|\tilde\Delta B^H_{\frac{i-1}{n}}\right| = 2^{H+\frac{1}{2}}T, \end{align*}

in probability. Taking into account of $H \gt \frac{1}{2}$ , we have

\begin{align*}\lim_{n \to \infty} V^n(X)_T = \frac{1}{2} \lim_{n \to \infty} U^n(X)_T = \alpha^2 \sigma^2 T, \end{align*}

and

\begin{align*}\lim_{n \to \infty} n^{H-\frac{1}{2}}\left(-2V^n(X)_T + U^n(X)_T\right) = 2\alpha\sigma^2 T\left(2^{H+\frac{1}{2}}-2\right), \end{align*}

in probability. Based on the above results, we can construct the consistent estimators for the parameters $\alpha$ and $\sigma$ :

(18) \begin{align} \hat\alpha = \frac{2 \left(2^{H+\frac{1}{2}}-2\right) V^n(X)_T}{n^{\hat H-\frac{1}{2}}\left(U^n(X)_T - 2V^n(X)_T\right)}, \\[-6pt] \nonumber\end{align}

(19) \begin{align} \hat\sigma = \sqrt{\frac{V^n(X)_T}{\hat\alpha^2 T + T n^{1-2\hat H} + 2 \hat\alpha T n^{\frac{1}{2} - \hat H}}}. \\[10pt] \nonumber \end{align}

We have thus far discussed the estimation of the parameters H, $\alpha$ , and $\sigma$ . In the remainder of this subsection, we provide two approaches for estimating the parameters m and $\theta$ .

5.1.1. Data-driven approach for estimating m and $\theta$

The parameters m and $\theta$ in the drift part can be estimated using the least squares method (see, for example, Hu et al., Reference Hu, Nualart and Zhou2019 and Neuenkirch and Tindel, Reference Neuenkirch and Tindel2014 for the idea). Namely, we would like to minimize the function

\begin{align*}L(m, \theta) = \sum_{i=1}^{[nT]} \left|X_{\frac{i}{n}} - X_{\frac{i-1}{n}} - \frac{1}{n}\left(m - \theta X_{\frac{i-1}{n}}\right)\right|^2. \end{align*}

Then, the least squares estimators $\hat m$ and $\hat \theta$ can be defined by $(\hat m, \hat \theta) = \mathrm{argmin} \ L(m, \theta)$ , and we can derive the explicit expressions for the estimators:

(20) \begin{equation} \hat m = n\frac{\left(\sum_{i=1}^{[nT]} \Delta X_{\frac{i-1}{n}} \right) \left(\sum_{i=1}^{[nT]} X^2_{\frac{i-1}{n}}\right) - \left(\sum_{i=1}^{[nT]} \Delta X_{\frac{i-1}{n}} \cdot X_{\frac{i-1}{n}} \right) \left(\sum_{i=1}^{[nT]} X_{\frac{i-1}{n}}\right)}{ [nT]\sum_{i=1}^{[nT]} X^2_{\frac{i-1}{n}} - \left(\sum_{i=1}^{[nT]} X_{\frac{i-1}{n}}\right)^2}, \end{equation}

(21) \begin{equation} \hat \theta = n\frac{\left(\sum_{i=1}^{[nT]} \Delta X_{\frac{i-1}{n}} \right) \left(\sum_{i=1}^{[nT]} X_{\frac{i-1}{n}} \right) - [nT]\sum_{i=1}^{[nT]}\Delta X_{\frac{i-1}{n}} \cdot X_{\frac{i-1}{n}}}{[nT]\sum_{i=1}^{[nT]} X^2_{\frac{i-1}{n}}- \left(\sum_{i=1}^{[nT]} X_{\frac{i-1}{n}}\right)^2}. \end{equation}

5.1.2. Judgment-driven approach for estimating m and $\theta$

As an alternative to the data-driven approach, the parameters m and $\theta$ can be estimated by ergodic-type estimators, which are consistent and do not require a high-frequency data structure, making them particularly straightforward to implement in practice (Hu et al., Reference Hu, Nualart and Zhou2019; Wang et al., Reference Wang, Xiao and Yu2023). Namely, the mean reversion rate $\theta$ can be estimated by

(22) \begin{equation} \hat\theta = \left(\frac{[nT]\sum_{i=1}^{[nT]} X_{\frac{i}{n}}^2 - \left(\sum_{i=1}^{[nT]} X_{\frac{i}{n}}\right)^2}{[nT]^2 \hat\sigma^2 \hat H \Gamma(2\hat H)}\right)^{-\frac{1}{2 \hat H}}.\end{equation}

The quantity $\frac{m}{\theta}$ is known as the long-term mean of $X_t$ , and it can be determined from an external resource (such as expert opinions or actuarial judgment on the long-term trend of $X_t$ ). Based on the estimation of the long-term mean, we can estimate the parameter m.

(23) \begin{equation} \hat m = \bar{X} \times \hat\theta, \end{equation}

where $\bar{X}$ denotes the actuarially judged long-term mean of $X_t$ . This approach is flexible enough to incorporate practitioners’ actuarial judgment on the long-term trend of the process.

We have presented two approaches above, and the choice between them depends on the nature of the data. For instance, interest rates are largely impacted by monetary policy, so the judgment-driven approach could be more appropriate. In contrast, mortality dynamics are less directly governed by policy interventions, making the data-driven approach more appropriate.

5.2. Model calibration for the bivariate process $(r_t, \mu_t)$

In this section, we apply the methods described in Section 5.1 to estimate the eleven parameters in the bivariate process $(r_t, \mu_t)$ given by Equation (4). The sequential algorithm is described in the following six steps.

Step 1: Apply the rescaled range analysis method to estimate the Hurst parameter $H_1$ for the interest rate process $r_t$ .

Step 2: Apply the corresponding equations in Section 5.1 to estimate the parameter $\alpha_1$ and the volatility parameter $\sigma_1$ for the interest rate process $r_t$ .

Step 3: Estimate the parameters $m_1$ and $\theta_1$ for the interest rate process $r_t$ based on ergodic-type estimators; i.e., Equations (22) and (23).

Step 4: Repeat Steps 1–2 to estimate the parameters $H_2$ , $\alpha_2$ , and $\sigma_2$ for the excess mortality process $\mu_t$ .

Step 5: Estimate the parameters $m_2$ and $\theta_2$ for the excess mortality process $\mu_t$ based on least squares estimators; i.e., Equations (20) and (21).

Step 6: Estimate the parameter $\rho$ using the following variance-adjusted estimator based on Lemma 3.1:

\begin{align*}\hat\rho = \frac{\sum_{i=1}^{[nT]} (e_i^{(r)} - \bar e^{(r)})(e_i^{(\mu)} - \bar e^{(\mu)}) \sqrt{\sigma_1^2 \alpha_1^2 n^{-1} + \sigma_1^2 n^{-2H_1}}\sqrt{\sigma_2^2 \alpha_2^2 n^{-1} + \sigma_2^2 n^{-2H_2}}}{\sqrt{\sum_{i=1}^{[nT]} (e_i^{(r)} - \bar e^{(r)})^2}\sqrt{\sum_{i=1}^{[nT]} (e_i^{(\mu)} - \bar e^{(\mu)})^2}\hat\sigma_1 \hat\sigma_2 \hat\alpha_1 \hat\alpha_2 n^{-1}}, \end{align*}

where the sequences $e_i^{(r)}$ and $e_i^{(\mu)}$ can be derived from

\begin{align*}e_i^{(r)} = r_{\frac{i}{n}} - r_{\frac{i-1}{n}} - \left(\hat m_1 - \hat\theta_1 r_{\frac{i-1}{n}}\right) \frac{1}{n}, \quad e_i^{(\mu)} = \mu_{\frac{i}{n}} - \mu_{\frac{i-1}{n}} - \left(\hat m_2 - \hat\theta_2 \mu_{\frac{i-1}{n}}\right) \frac{1}{n}, \end{align*}

with $\bar e^{(r)} = \frac{1}{[nT]} \sum_{i=1}^{[nT]} e_i^{(r)}, \ \bar e^{(\mu)} = \frac{1}{[nT]} \sum_{i=1}^{[nT]} e_i^{(\mu)}$ .

5.3. Calibration of pricing parameters

The previous section addressed the parameter estimation problem for the excess mortality and interest rate models under the physical measure. In this section, we present a general calibration procedure for estimating the pricing parameters under a selected risk-neutral pricing measure. The proposed methodology will be implemented later using market data from the Swiss Re Vita VI bond.

We assume that the mortality bond is a par value bond, which means $P_{0}^{\mathrm{MB}}(c) = F$ and the risk premium is reflected in the coupon rate c. In other words, the coupon rate c should be the risk-free rate at bond issuance plus a premium spread such that the bond price calculated from Equation (17) equals F. Setting $P_{0}^{\mathrm{MB}} = F$ , we can derive the formula of c from Equation (17) to obtain

(24) \begin{equation} c = \frac{1 - P(0, T) + \mathbb{E}^{\mathbb{Q}}\left[\mathrm{PRF}(\tilde \mu) \times e^{-\int_0^T r_s {\mathop{}\! {\mathop{}\!\mathrm{d}}s}} \big|\mathcal{F}_0^r\right]}{\frac{1}{m}\sum_{0 \lt t_k \leq T}P(0, t_k)}.\end{equation}

The MLS traded in the market generally provides information about the coupon rate $c_{\mathrm{obs}}$ , expected loss $EL_{\mathrm{obs}}$ , probability of attachment $PA_{\mathrm{obs}}$ , and premium spread $PS_{\mathrm{obs}}$ . The general principle is to match the mathematical equations derived from the stochastic models with the market information. The calibration process for the MLS involves numerical solutions for the following parameters: attachment point a, detachment point b, risk premiums $\gamma_1, \gamma_2, \eta_1, \eta_2$ . In general, we can follow the algorithm below to estimate these parameters.

Step 1: Obtain simulated mortality rate trajectories over the MLS term from the mortality model under physical measure $\mathbb{P}$ . For our proposed model, we can simulate the excess mortality rates from the model (4) and recover the mortality rates by adding back the selected baseline mortality rates.

Step 2: Calculate the mortality index $\tilde\mu$ based on the MLS design as discussed in Section 4.2 for all simulated trajectories obtained in Step 1 and derive the estimated distribution of $\tilde\mu$ . If the probability of attachment information is available on the market, the attachment point a can be estimated by the $100(1-PA_{\mathrm{obs}})$ th percentile of the distribution of $\tilde\mu$ . If the expected loss information is available, we can obtain an estimation of the detachment point b through a selected optimization method such that

\begin{align*}EL_{\mathrm{obs}} = \mathbb{E}[\mathrm{PRF}(\tilde \mu)]. \end{align*}

Step 3: For candidate values of the risk premium parameters, simulations are performed for the interest rate and the excess mortality rate under the pricing measure $\mathbb{Q}$ so that the right-hand side of Equation (24) can be numerically approximated. Then, an optimization approach can be implemented to minimize the absolute differences between the observed coupon rates and the approximated coupon rates from numerical simulations. The optimization process is similar if the premium spread information is available, as it is essentially the difference between the coupon rate and the risk-free interest rate.

The above general methodology can be applied for many types of MLS, but the specific formulas and selection of the optimization procedure will be customized and provided for a real-world mortality bond in Section 6.2.2.

6.1. Data description

The mortality and interest rate data used in our numerical illustration are obtained from the HMDFootnote ² and the FRED,Footnote ³ respectively. The two series are described in detail below:

• • Mortality data: Weekly all-age unisex death counts $D_t$ and the corresponding exposure estimates $E_t$ are taken directly from the short-term mortality fluctuations (STMF) series of the HMD. The excess mortality rate in week t is defined as the deviation of the observed weekly death rate from a pre-specified baseline: $\mu_t = \frac{D_t}{E_t} - \mu_t^{(b)}$ , where $\mu_t^{(b)}$ denotes the baseline death rate. To account for seasonal variations in weekly mortality, the baseline for the m-th week of the year is estimated as the average of the observed weekly death rates for that same calendar week across the years 2015–2019; i.e., the five pre-pandemic years in our dataset.

• • Interest rate data: Weekly interest rates are taken from the WTB3MS series of the FRED, which reports the 3-month U.S. Treasury bill rate on a secondary-market discount basis. As this series is directly observable, no further construction is required other than expressing the rates as decimal values $r_t$ rather than percentages.

The sample period for both series spans January 2015 to December 2024, comprising a total of 521 weekly observations. This window is chosen to align with the availability of U.S. weekly mortality data in the STMF series, and the same period is applied to the interest rate data to ensure consistency.

Table 1.

Estimated parameters of the proposed model calibrated to weekly mortality and interest rate data from 2015 to 2024.

Table 1. Long description

A table with two rows and seven columns comparing estimated parameters of a model calibrated to weekly mortality and interest rate data from 2015 to 2024. The columns are labeled H, alpha, sigma, m, theta, and rho. The row labels are ri and mu i. Row 1: H, 0.85723; alpha, 0.24925; sigma, 0.00918; m, 0.01584; theta, 0.37897; rho, -0.47591. Row 2: H, 0.78373; alpha, 0.32554; sigma, 0.00213; m, 0.00072; theta, 1.18381.

6.2. Calibration results

This section presents the calibration results for our proposed model (4) and its risk neutral version (9). While a general overview of the calibration procedure for the risk premium parameters is provided in Section 5.3, this section offers a more detailed explanation of the procedure customized for the data described in Section 6.1.

6.2.2. Calibration under measure $\mathbb{Q}$ for pricing parameters

We use real market data to calibrate the pricing parameters under the risk-neutral measure $\mathbb{Q}$ , as described in Theorem 3.4. In particular, we calibrate the pricing parameters to match the price of the Vita Capital VI bond, issued in July 2021. The relevant pricing information of this bond, collected from Artemis,Footnote ⁵ is reported in Table 2. Based on the available market information, we assume that the Vita VI bond is a 5-year par value bond that makes annual payments following Equation (16) at a specified coupon rate of 3%. The mortality index that determines the bond’s principal reduction factor is the average mortality rate observed over each year during the 5-year term.

Table 2.

Pricing information of the Vita Capital VI bond.

Table 2. Long description

A table displaying pricing information of the Vita Capital VI bond. The table has three columns and one row. The columns are labeled Term, Probability of first loss, Expected loss, and Coupon rate. The row provides data for a 5-year term: Probability of first loss is 1.06 percent, Expected loss is 0.75 percent, and Coupon rate is 3 percent.

The general methodology in Section 5.3 can be applied to find other unknown parameters of this specific mortality bond. Step 1 and Step 2 in Section 5.3 can be implemented to determine the attachment point and the detachment point. To match the coupon rate, the evaluation of Equation (17) is not straightforward, and it may not be easy to find a closed formula due to the challenge of deriving the conditional expectation of a functional that involves extremal functionals with respect to fBm. However, we will perform the simulation-based optimization procedure as described in Step 3 of Section 5.3.

As we have limited data information, to reduce the dimensionality of the parameter’s space for optimization, we mandate that the risk premiums for the fBm part are zero, that is, $\eta_1 = \eta_2 = 0$ . To determine the risk premium $\gamma_1$ , we can implement an optimization procedure to match the price of the zero coupon bond under the measure $\mathbb{Q}$ with the one determined by the target interest rate i. That is, the optimal value of the risk premium $\gamma_1$ is given by: $\mathrm{argmin}_{\gamma_1} |P(0, T)_{m_1'} - (1+i)^{-5}|,$ where $P(0, T)_{m_1'}$ is the price of the zero coupon bond from Theorem 4.2 under the probability measure $\mathbb{Q}$ with the parameter $m_1' = m_1 + \alpha_1 \sigma_1 \gamma_1 + \sigma_1 \eta_1$ given in Theorem 3.4. The interest rate i measures the expected yield rate of a 5-year corporate bond for pricing purposes. We choose $i=2.57\%$ , which is the monthly average AAA bond (annualized) spread when the mortality bond was issued. We adopt the fmincon optimization method in Matlab to obtain the optimal value of the parameter $\gamma_1$ as $\hat\gamma_1 = 3.8701$ .

The attachment point can be found numerically by finding a value a such that the probability of attachment (i.e., the probability of first loss (PFL)) is 1.06%: $P\left(\mathrm{PRF}(\tilde{\mu}) \gt 0\right) = 1.06\%.$ Note that $\mathrm{PRF}(\tilde{\mu}) \gt 0$ if and only if $\{\tilde{\mu}_{t_k} \gt a | k = 1, \ldots, 5\} \neq \emptyset$ . We can simulate 10,000 trajectories of excess mortality according to the SDE under the pricing measure $\mathbb{Q}$ and add the baseline mortality rate (i.e., the expected death rate $\mathbb{E}[D_t/E_t]$ ) to the simulated excess mortality. In this way, we obtain the 10,000 simulated weekly mortality rates $\mu_t$ . The principal reduction factor (PRF) is then determined as

\begin{align*}\mathrm{PRF}(\tilde{\mu}) = \min\left(1, \sum_{t_k= 1}^5 \frac{(\tilde{\mu}_{t_k} - a)_+ - (\tilde{\mu}_{t_k} - b)_+}{b - a} \right), \end{align*}

where $\tilde{\mu}_{t_k}$ is the average mortality rate for year $t_k$ , determined by the simulated weekly mortality rates and computed as $\tilde{\mu}_{t_k} =\frac{\sum_{m=1}^{52}\mu_{52\times t_k +m}}{52}$ . The $98.94$ th percentile of all values $\{\tilde{\mu}_{t_k}\;:\; k=1,\ldots,5\}$ will be the estimated value of the attachment point $\hat a$ .

Given the expected loss data information, we can set up the equation $\mathbb{E}\left[\mathrm{PRF} (\tilde{\mu})\Big|\mathrm{PRF}(\tilde{\mu}) \gt 0 \right] \times \Pr\left(\mathrm{PRF}(\tilde{\mu}) \gt 0\right) = 0.75\%.$ The detachment point b can be estimated by solving the above equation via minimizing the function

\begin{align*}f(b) = \left|\frac{1}{N}\sum_{i=1}^N \min\left(1, \sum_{t_k= 1}^T \frac{(\tilde{\mu}_{i, t_k} - a)_+ - (\tilde{\mu}_{i, t_k} - b)_+}{b - a} \right) - \frac{0.75\%}{1.06\%}\right|, \end{align*}

where $\tilde{\mu}_{i, t_k}$ is the ith simulated value of $\tilde{\mu}_{t_k}$ that satisfies $\tilde{\mu}_{i, t_k} \gt a$ , and N is the number of simulated trajectories such that the attachment point is triggered.

Our last step is to find the value of the risk premium $\gamma_2$ such that the coupon rate 3% satisfies the pricing Equation (24). To achieve this purpose, we implement the grid search method for the target value $\gamma_2$ in the interval [0, 2] with step size $0.001$ , such that

\begin{align*}3\% = \frac{1 - P(0, T) + \frac{1}{N} \sum_{i=1}^N \min\left(1, \sum_{t_k= 1}^T \frac{(\tilde{\mu}_{i, t_k} - a)_+ - (\tilde{\mu}_{i, t_k} - b)_+}{b - a} \right) \times e^{- \sum_{t_k=1}^5 \sum_{m=1}^{52} r_{i, 52t_k + m}}}{\frac{1}{m}\sum_{0 \lt t_k \leq T}P(0, t_k)} \end{align*}

holds numerically, where $r_{i, 52t_k + m}$ are the simulated short rates under pricing measure $\mathbb{Q}$ . The optimal value for the risk premium $\gamma_2$ is finally found to be $\hat\gamma_2 = 1.0620$ .

6.3. Baseline analysis

After calibrating the model and pricing parameters to real mortality data and market information, we now analyze a 5-year catastrophe mortality bond structured similarly to the Vita VI bond. We set the face value to $F = 100$ and assume that the bond makes annual payments according to the formula (16). The mortality index that determines the bond’s principal reduction factor is defined as the yearly average mortality rate over the 5-year term. We use simulated paths of mortality and interest rates under measure $\mathbb{Q}$ to determine the coupon rate, and under measure $\mathbb{P}$ to compute the present value of the bond’s payouts at issuance (i.e., time 0).

Figure 3.

Simulated paths of weekly mortality rates (left) and interest rates (right) with 95% confidence intervals under $\mathbb{Q}$ -measure.

6.3.2. Principal repayment

We now investigate the payments of the catastrophe mortality bond using our calibrated model under the physical measure $\mathbb{P}$ . Specifically, we analyze the distribution of principal repayment and total payouts based on simulated mortality and interest rates. Figure 4(a) presents the distribution of simulated principal repayment values at time 0, computed under the physical measure $\mathbb{P}$ . Each value represents the time-0 present value of the principal repayment paid at maturity, discounted using the corresponding simulated interest rate paths. The repayment amount is determined by the principal reduction factor, which is calculated using the simulated mortality rate paths over the bond’s 5-year term.

Table 3.

Risk measures for the present value (PV) of principal repayment and the present value of total bond payouts under the physical measure $\mathbb{P}$ .

Table 3. Long description

A table comparing risk measures for the present value of principal repayment and total bond payouts under the physical measure. The table has two rows and seven columns. The columns are labeled Mean, Std, VaR 5 percent, CTE 5 percent, VaR 1 percent, and CTE 1 percent. The first row represents the present value of principal repayment with values: Mean 80.34, Std 7.91, VaR 5 percent 73.43, CTE 5 percent 61.48, VaR 1 percent 67.52, and CTE 1 percent 19.93. The second row represents the present value of total bond payouts with values: Mean 105.37, Std 8.30, VaR 5 percent 97.36, CTE 5 percent 85.40, VaR 1 percent 90.74, and CTE 1 percent 44.76.

Figure 4.

Histogram of simulated present values of principal repayment (left) and total payouts (right) at time 0, computed under the physical measure $\mathbb{P}$ using simulated mortality and interest rate paths.

The main histogram of Figure 4(a) shows the distribution of discounted principal repayments across all simulated mortality and interest rate scenarios. A large proportion of the mass is concentrated near full repayment, reflecting typical cases where the mortality index remains below the attachment point and investors recover the entire principal. To better examine tail risks, the figure includes a zoom-in on the lower repayment region, specifically values between 10 and 60. This region corresponds to extreme mortality scenarios, where high mortality rates trigger significant reductions in principal repayment. Lastly, a small but noticeable proportion of the mass is concentrated at zero repayment, indicating the likelihood of the most severe cases.

Table 3 summarizes key numerical risk measures of the simulated principal repayment values at time 0. Specifically, we report the mean, standard deviation (Std), Value-at-Risk at the 5% level (VaR $_{5\%}$ ) and 1% level (VaR $_{1\%}$ ), and Conditional Tail Expectation at the 5% level (CTE $_{5\%}$ ) and 1% level (CTE $_{1\%}$ ). The mean repayment amount is approximately $80.34$ , indicating that investors typically recover a substantial portion of the principal under common mortality scenarios. The tail risk measures offer a deeper insight into extreme scenarios. The VaR $_{5\%}$ value of $73.43$ indicates that there is a 5% probability that principal repayment will fall below this threshold, while the CTE measures emphasize the severity of these tail scenarios. The CTE $_{5\%}$ at $61.48$ and the substantially lower CTE $_{1\%}$ at $19.93$ clearly show that extreme mortality events can lead to substantial losses.

6.4. Sensitivity analysis

In this subsection, we conduct a series of sensitivity analyses to better understand how the bond’s outcomes respond to changes in input assumptions. We examine key quantities, including loss metrics, payout measures, and the fair coupon rate. Three categories of parameters are considered: (i) bond design features, such as attachment and detachment points; (ii) model parameters estimated under the physical measure; and (iii) pricing parameters calibrated under the risk-neutral pricing measure. Each scenario builds on the baseline results presented earlier. Unless otherwise specified, all other inputs and bond settings remain fixed at their baseline values.

6.4.1. Sensitivity to bond parameters

We first analyze how bond structural parameters, namely, the attachment point, detachment point, and bond duration, influence the fair coupon rate of a catastrophe mortality bond.

Figure 5 presents a heatmap of the fair coupon rate as a function of the attachment and detachment points, expressed as multiples of the mortality experience in the bond issue year 2024. As expected, the coupon rate decreases monotonically as either the attachment point or the detachment point increases. A higher attachment point reduces the likelihood of triggering principal reduction, while a higher detachment point limits the severity of potential principal losses.

More interestingly, the sensitivity shown in Figure 5 is not uniform across the parameter space. We observe that changes in the detachment point have a more pronounced effect when the attachment point is low. As the attachment point increases, the fair coupon rate becomes less sensitive to further increases in the detachment point. This pattern reflects a diminishing marginal effect of expanding the reduction coverage once the probability of triggering any principal reduction becomes sufficiently small.

Table 4 reports pricing metrics for catastrophe mortality bonds with maturity terms from 1 to 5 years, including the principal repayment’s PFL, conditional expected loss (CEL), expected loss (EL), and the fair coupon rate. All values are based on 10,000 simulations under the calibrated risk-neutral measure $\mathbb{Q}$ . As the bond term shortens, both the probability and severity of losses decline, reflecting reduced exposure to extreme mortality events. This further leads to lower risk premiums, as reflected in steadily decreasing coupon rates.

Table 4.

Loss metrics of the principal repayment and fair coupon rates for catastrophe mortality bonds with varying maturity terms, calculated under the risk-neutral measure $\mathbb{Q}$ .

Table 4. Long description

The table presents pricing metrics for catastrophe mortality bonds with maturity terms ranging from 1 to 5 years. It includes columns for principal repayments PFL, conditional expected loss (CEL), expected loss (EL), and the fair coupon rate. The table has five rows, each corresponding to a different maturity term: 5 years (baseline), 4 years, 3 years, 2 years, and 1 year. Row 1: 5 years (baseline), PFL 1.06 percent, CEL 70.75 percent, EL 0.75 percent, Coupon rate 5.68 percent. Row 2: 4 years, PFL 0.79 percent, CEL 70.74 percent, EL 0.56 percent, Coupon rate 5.55 percent. Row 3: 3 years, PFL 0.36 percent, CEL 68.56 percent, EL 0.25 percent, Coupon rate 5.36 percent. Row 4: 2 years, PFL 0.15 percent, CEL 58.20 percent, EL 0.08 percent, Coupon rate 5.13 percent. Row 5: 1 year, PFL 0.00 percent, CEL N/A, EL N/A, Coupon rate 4.83 percent.

Figure 5.

Heatmap showing the fair coupon rate as a function of the attachment point and detachment point, calculated under the risk-neutral pricing measure $\mathbb{Q}$ .

The fair coupon drops from 5.68% for a 5-year bond to 4.83% for a 1-year bond. For the 1-year bond, no loss events were observed, so CEL and EL are not reported. The corresponding coupon rate thus reflects only interest rate risk. Although the CEL remains nearly unchanged between the 4-year and 5-year bonds, the expected loss decreases slightly, consistent with the reduction in PFL from 1.06% to 0.79%. These results showcase how bond duration affects loss likelihood, severity and coupon rate of catastrophe mortality bonds.

6.4.2. Sensitivity to model parameters

We now investigate the sensitivity of bond outcomes to various model parameters. In particular, we consider the following scenarios:

• • Scenario 1: The instantaneous correlation parameter is set to zero.

• • Scenario 2: The Hurst parameters for both $r_t$ and $\mu_t$ are set to 0.5.

• • Scenario 3: The variance parameter of the interest rate process is doubled.

• • Scenario 4: The variance parameter of the excess mortality rate process is doubled.

• • Scenario 5: The pricing parameter associated with $r_t$ is increased by 50%.

• • Scenario 6: The pricing parameter associated with $\mu_t$ is increased by 50%.

Table 5 reports the fair coupon rate of the bond under each scenario, along with summary statistics and loss metrics for the present value of total bond payouts. For reference, the corresponding quantities under the baseline scenario are also included. As in the baseline case, the fair coupon rates are determined using simulations under the risk-neutral measure $\mathbb{Q}$ , while the payout risk measures are computed under the physical measure $\mathbb{P}$ .

Table 5.

Sensitivity analysis results for the scenarios under consideration. Each row reports summary statistics of the present value of total bond payouts (mean, standard deviation, VaR, CTE), loss metrics of the principal repayment (PFL, CEL, EL), and the fair coupon rate.

Table 5. Long description

The table presents sensitivity analysis results for various scenarios, focusing on the present value of total bond payouts and associated metrics. It includes seven rows and nine columns. The columns are labeled as Scenario, Mean, Std, VaR, CTE, PFL (%), CEL (%), EL (%), and Coupon rate (%). The row labels include Baseline, Scenario 1, Scenario 2, Scenario 3, Scenario 4, Scenario 5, and Scenario 6. Each row provides values for the mean, standard deviation, VaR, CTE, PFL, CEL, EL, and coupon rate for each scenario. Notable trends include variations in the mean, standard deviation, and coupon rates across different scenarios, with Scenario 4 showing the highest standard deviation and Scenario 5 the highest mean. The table provides a detailed comparison of how different scenarios affect the bond payouts and associated metrics.

Scenarios 1 and 2 examine the effects of removing key dependence structures between mortality and interest rates. When correlation is removed in Scenario 1, the coupon rate increases noticeably and payout risk measures shift relative to the baseline, indicating that short-term co-movement between the two processes materially influences pricing. In contrast, removing LRD in Scenario 2 results in a small reduction in the coupon rate and noticeably lower tail risk, suggesting that long-memory effects remain important for characterizing extreme outcomes. These findings confirm the need of incorporating both correlation and long-memory features into the joint modeling framework, as each plays a distinct role in shaping risk and return in catastrophe mortality bonds.

Scenarios 3 and 4 examine the effects of higher volatility in the interest rate process and the excess mortality rate process, respectively. Doubling the variance of interest rate in Scenario 3 leads to a noticeable increase in the coupon rate and greater variability in the bond payouts, indicating that interest rate volatility has a meaningful influence on bond pricing and return uncertainty. In contrast, doubling the variance of mortality rate in Scenario 4 results in a dramatic increase in tail risk, accompanied by a significantly higher coupon rate. This result suggests that mortality volatility is a comparatively stronger driver of both pricing and risk in catastrophe mortality bonds.

Scenarios 5 and 6 assess the impact of higher risk premiums for interest rate risk and mortality risk, respectively. In Scenario 5, increasing the interest rate risk premium leads to a noticeable increase in the coupon rate, reflecting the higher compensation required by investors for bearing additional interest rate risk. In Scenario 6, increasing the mortality risk premium by the same proportion produces a more moderate increase in the coupon rate. Overall, interest rate risk premiums have a stronger influence on investor compensation than mortality risk premiums. Lastly, the loss metrics under both scenarios remain unchanged relative to the baseline because these metrics are irrelevant to how risks are compensated under the risk-neutral measure $\mathbb{Q}$ .