2.1. Stochastic mortality modeling with Volterra model

In this subsection, we introduce basic definitions and concepts of stochastic mortality and the Volterra model. First, we consider a probability space $(\Omega, \mathcal{F}, \mathbb{F}, \mathbb{P})$ , where the filtration $\mathbb{F} = (\mathcal{F}_t)_{t \geq 0}$ satisfies the usual conditions and represents the information flow related to the mortality process. And $\mathbb{P}$ is the usual physical probability measure under which the mortality model is defined. Let N(s) denote an $\mathbb{F}$ -adapted counting process with associated stopping time $\tau_x$ , representing the future lifetime of a policyholder age x at time $s=0$ (with t denote the time of entrance). The mortality intensity process $\mu_x(s)$ , governing the counting process evolution, resides in a sub-filtration $\mathbb{G} \subseteq \mathbb{F}$ that excludes individual death time information. $\mathbb{G}$ may be interpreted as representing the evolution of mortality intensity $\mu$ , but not the actual mortality experience of an individual. Given a realization of state $\omega \in \Omega$ and historical path for intensity up to time t. For a homogeneous cohort of age x, the conditional survival probability satisfies:

\[P(N(s+\Delta s)-N(s)\gt 0|\mathcal{F}_s) \approx \mu_x(s)\Delta s \quad \text{for} \quad \tau_x \gt s. \]

This stochastic formulation fundamentally differs from deterministic models through its explicit incorporation of intensity path dependence. Conditional on trajectory realization over [s, T], the counting process follows an heterogeneous Poisson law:

\[P(N(T)-N(s)=k|\mathcal{F}_s \vee \mathcal{G}_T) = \frac{\left(\int_s^T \mu_x(r) dr\right)^k \exp\left\{-\int_s^T \mu_x(r) dr\right\}}{k!}, \]

and setting $k=0$ gives the usual deterministic survival function.

Then we refine the baseline Gompertz mortality model by incorporating LRD to better reflect real-life mortality patterns. Recognizing that an individual’s mortality rate is influenced not only by current factors but also by historical experiences and environmental influences, we propose a dynamic approach to model mortality rates. The Volterra mortality model is first introduced and extensively studied in Wang et al. (Reference Wang, Chiu and Wong2021), based on key results and theorems from Biffis (Reference Biffis2005) and Abi Jaber et al. (Reference Abi Jaber, Larsson and Pulido2019b). Following this framework, we express the mortality rate as

(2.1) \begin{equation} \mu_x(s) = M_x(s, Y_s) = m_x(s) + \eta Y_{x}(s), \end{equation}

where $m_x(s)$ represents the baseline mortality rate, which increases over time and with the age of the individual. The baseline mortality rate is modeled using the Gompertz function as follows:

(2.2) \begin{equation} m_{x}(s) = \frac{1}{b} e^{\frac{x+s-m}{b}}, \end{equation}

where b is the scaling factor for the Gompertz function, determining the rate at which mortality increases, $m\gt 0$ denotes the modal age at death. The process $\{Y_x(s)\}_{s \geq 0}$ follows a stochastic Volterra integral equation:

(2.3) \begin{equation} Y_x(s) = Y_x(0) + \int_0^s K(s-r)b(Y_x(r)) \, dr + \int_0^s K(s-r)\sigma(Y_x(r)) \, dW(r), \end{equation}

where $Y_x(0) $ is the initial value of the process, $ K(s - r) $ is a kernel function that determines the weight of past events at time r on the current time s. This kernel introduces memory effects, enabling the model to capture LRD by weighting historical influences on the current state. And W(s) is a one-dimensional standard Brownian motion, which introduces stochasticity into the process. Here, b(y) and $ \sigma(y) $ are defined as follows:

\[ b(y) = b_0 + b_1 y, \quad \sigma(y) = \sigma_0 + \sigma_1 y, \]

where $ b_0 $ , $ b_1 $ , $ \sigma_0 $ , and $ \sigma_1 \in \mathbb{R} $ are constants. These equations define a dynamic mortality model that incorporates LRD, marking an extension to the traditional tontine model.

The kernel function $K(s-r)$ encodes historical memory. It satisfies:

\[ K \in L^2_{\text{loc}}(\mathbb{R}_+,\mathbb{R}), \quad \int_0^h K(t)^2\,dt = O(h^\gamma), \quad \int_0^T (K(t+h)-K(t))^2 \,dt = O(h^\gamma), \]

for some $\gamma \in (0,2]$ and any $T\lt \infty$ . This ensures well-posedness of the Volterra process. If K is constant, $(Y_s)_{s\ge0}$ reduces to a standard affine SDE. By contrast, a non-constant kernel enables LRD. A key example is the fractional kernel:

\[ K(t) = \frac{c\,t^{\alpha-1}}{\Gamma(\alpha)}, \]

where $c\gt 0$ and $\Gamma(\cdot)$ is the Gamma function, and where $\alpha = H + 1/2$ and $H \in [0.5,1]$ is the Hurst parameter H. This kernel assigns decaying weights to past mortality shocks via a power-law structure, aligning with empirical observations of persistent memory in mortality data.

Chen et al. (Reference Chen, Hieber and Klein2019) instead introduce a shock variable $\varepsilon$ for systematic risks, uniformly shifting mortality rates across all ages. While effective for simplicity, this approach neglects age-specific deviations and path-dependent memory. However in our case, since deviations from the baseline mortality function (or superimposed random noise) are modeled directly through a Volterra process underlying our model and LRD is implied when a general kernel (non-constant) K is adopted, we do not consider any forms of shocks separately. Once K is non-constant, LRD arises naturally. In this way, the model extends classical tontine or affine mortality frameworks. The choice of modeling the deviations in mortality rates using the stochastic Volterra process $Y_x(s)$ is motivated by the need to capture both the stochastic fluctuations and the LRD in mortality dynamics. The kernel function $K(s-r)$ allows the process to incorporate historical dependencies, meaning that the evolution of $Y_x(s)$ depends on its entire past trajectory. This is essential for accurately modeling mortality rates, as past mortality experiences and environmental factors can have a persistent impact on future mortality.

2.2. Tontine structure under Volterra mortality

In this subsection, we formalizes the retirement income tontine structure under the Volterra mortality model introduced in Subsection 2.1. We first define the tontine payout mechanism for a homogeneous cohort of policyholders, where mortality dynamics are governed by the stochastic intensity $\mu_x(s)$ in (2.1) and its Volterra-driven deviation process $Y_x(s)$ in (2.3). The framework is then extended to multiple cohorts, each with independent mortality trajectories. All formulations are conditioned on the filtered probability space $(\Omega, \mathcal{F}, \mathbb{F}, \mathbb{P})$ , where $\mathbb{F}$ captures mortality experience and financial market information.

For a homogeneous cohort of n policyholders, let t denote the time of entrance and aged x at entry, the mortality intensity process follows the Volterra model as defined in formula (2.1). Consider the indicator random variable $\mathbb{I}_{\{\tau_{x} \gt s\}}$ , the force of mortality of each member is $\mu_{x}(s)$ in this cohort. Furthermore, we introduce the notation:

\[ S_x(s) = \mathbb{E}\left[ \mathbb{I}_{\{\tau_x \gt s\}} \big| \mathcal{G}_{t+s} \right] = \exp\left( -\int_0^{s} \mu_x(r)dr \right), \]

where $ \mathcal{G}_{t+s} $ is the sigma-algebra generated by $ \{ \mu_{x}(r) \mid r \leq s \} $ , which contains all the information regarding the systematic mortality risk up to time s starting from age x. And the number of surviving participants N(s) follows a doubly stochastic binomial distribution:

\[ N(s) \big| \mathcal{G}_{t+s} \sim \text{Binomial}\left( n, \ S_x(s) \right).\]

Following the tontine structure and optimal payout problems that was initially proposed in Milevsky and Salisbury (Reference Milevsky and Salisbury2015), we adopt the similar retirement income scheme. The continuous benefit stream for the j-th surviving policyholder in the homogeneous cohort is:

(2.4) \begin{equation} b_j(s) = \frac{n w_j d(s)}{N(s)} \mathbb{I}_{\{\tau_x \gt s\}}, \end{equation}

where $j = 1, 2, \ldots, n$ and $w_j$ denotes initial contribution and d(s) is the payout rate per unit wealth which is uniform among all policyholders over the homogeneous cohort.

Next we extend the framework to L heterogeneous cohorts with distinct initial wealth levels and entry ages, where $n_i$ denotes the number of participants in cohort $i \in {1,\ldots,L}$ , yielding total initial pool size $n = \sum_{i=1}^L n_i$ . Members within each cohort are homogeneous in age and wealth, characterized by entry age $x_i$ (with t denote the time of entrance) and initial wealth $w_i$ per member in cohort i. Each cohort’s mortality intensity $\mu_{x_i}(s)$ follows:

(2.5) \begin{equation} \mu_{x_i}(s) = m_{x_i}(s) + \eta Y_{x_i}(s), \end{equation}

where $m_{x_i}(s)$ represents the Gompertz mortality and $Y_{x_i}(t)$ solves the independent Volterra equation:

(2.6) \begin{equation} Y_{x_i}(s) = Y_{x_i}(0) + \int_0^s K(s-r) b(Y_{x_i}(r)) dr + \int_0^s K(s-r) \sigma(Y_{x_i}(r)) dW_i(r), \end{equation}

where we further assume that $(W_1(t), W_2(t),\ldots,W_L(t))^{T}$ defines n-independent Brownian motions over the complete probability space.

Then the total survivors are $N(s) = \sum_{i=1}^L N_i(s)$ , where $N_i(s)$ counts survivors in cohort i. From the perspective of any surviving policyholder in the ith cohort, since $\mathbb{I} _{\{\tau_{x_i}\gt s\}}$ , the number of pool members $N_i(s)$ is at least 1, and the number of survivors N(s) conditional on the joint filtration $\mathcal{G}_{t+s}^{(1)} \vee \cdots \vee \mathcal{G}_{t+s}^{(L)}$ follows a Poisson-binomial distribution with probability vector

\[ \underbrace{S_{x_1}(s),\ldots,S_{x_1}(s)}_{n_1}, \ldots, \underbrace{S_{x_i}(s),\ldots,S_{x_i}(t)}_{n_i-1}, \ldots, \underbrace{S_{x_L}(s),\ldots,S_{x_L}(s)}_{n_L } , \]

of length $\sum_{i=1}^{L} n_i-1$ . Here, $S_{x_i}(s)$ is the conditional survival probability for cohort i:

\[ S_{x_i}(s) = \mathbb{E}\left[ \mathbb{I}_{\{\tau_{x_i} \gt s\}} \big| \mathcal{G}_{t+s}^{(i)} \right]=\exp\left( -\int_0^s \mu_{x_i}(r)\,dr \right), \]

where $\mathcal{G}_{t+s}^{(i)}$ is the sigma-algebra generated by $\{\mu_{x_i}(r)| r \leq s\}$ . Assuming independence between cohort-specific filtrations, N(s) reduces to the sum of L independent binomial random variables.

Following Milevsky and Salisbury (Reference Milevsky and Salisbury2016), let $d_i(s)$ denote the withdrawal rate for cohort i, reflecting the unique financial strategies or risk profiles adopted by each cohort, where $i=1,2,\ldots,L$ . The benefit adjusts for heterogeneity through tuning parameters $\pi_i$ :

(2.7) \begin{equation} b_{x_i}(s) \,:\!=\, \frac{\pi_i w d_i(s) }{N(s)} \mathbb{I}_{\{\tau_{x_i} \gt s\}},\end{equation}

where $w=\sum_{i=1}^{L}n_i w_i=\sum_{j=1}^{n} w_j$ represents the total initial wealth pooled from all participants, forming the basis of the collective fund. This variable can flexibly manage the allocation of funds to accommodate the specific needs and preferences of different groups. In addition, $\pi_i$ signifies the participation rate of individual in cohort i, which determines the individual’s proportional share within the total pooled wealth. Together, these parts form the basis for a comprehensive and equitable risk-sharing arrangement.

2.3. Actuarial fairness and f-value

In risk-sharing mechanisms, fairness is important for ensuring that all participants can benefit reasonably. Here, we formalize the concepts of individual actuarial fairness and collective actuarial fairness within the tontine model. First, we draw on the research of Milevsky and Salisbury (Reference Milevsky and Salisbury2015) and Chen et al. (Reference Chen, Hieber and Rach2021) regarding the individual fairness of all participants within a tontine and mutual aid domains. Subsequently, we investigate collective fairness, building on the work of Chen and Rach (Reference Chen and Rach2023) regarding collective fairness issues.

In an idealized model of expected actuarial equilibrium, the premiums paid by each member should adequately cover the expected amounts payable to other members upon their death. We fix a probability space $ (\Omega, \mathcal{F}, \mathbb{F}, \mathbb{P}) $ , where the filtration $\mathbb{F} = (\mathcal{F}_t)_{t \geq 0}$ satisfies all the usual conditions and consists of flow of information that includes the mortality experience. Let r be the risk free rate. Based on formula (2.7), for each policyholder in L distinct groups, $i = 1, 2,\ldots,L$ , we calculate the expected present value (EPV) of their periodic payments as follows:

(2.8) \begin{align}P_0^i &= \int_0^\infty e^{-rs}\mathbb{E} \left[ \frac{\pi_i w d_i(s)}{N(s)} \mathbb{I}_{\{\tau_{x_i} \gt s\}} \middle | \mathcal{F}_t\right] ds \notag \\[6pt] &= \int_0^\infty e^{-rs} \mathbb{E} \left[ \frac{\left( \sum_{l=1}^L n_l w_l \right) \pi_i d_i(s)}{\sum_{l=1}^L N_l(s)} \mathbb{I}_{\{\tau_{x_i} \gt s\}} \middle | \mathcal{F}_t \right] ds \notag \\[6pt] &= \int_0^\infty e^{-rs} \mathbb{E} \left[ S_{x_i}(s) \mathbb{E} \left[ \frac{\left( \sum_{l=1}^L n_l w_l \right) \pi_i d_i(s)}{\sum_{l=1}^L N_l(s) } \middle | \tau_{x_i}\gt s\right] \middle | \mathcal{F}_{t} \right] ds \notag \\[6pt] &= \int_0^\infty e^{-rs} \mathbb{E} \left[ S_{x_i}(s) \mathbb{E} \left[ \frac{\left( \sum_{l=1}^L n_l w_l \right) \pi_i d_i(s)}{\sum_{l=1}^L N_l(s) } \middle | \mathcal{G}_{t+s}\right] \middle | \mathcal{F}_{t} \right] ds.\end{align}

Definition 2.1. A tontine with homogeneous members is deemed fair if each member’s initial wealth $w_j$ can be expressed as

\[ w_j = P_0^{j} \quad \text{for } j=1,2,\ldots,n. \]

Individual fairness concentrates on balancing the benefits and contributions of each participant, ensuring each member’s contributions align with their received earnings. In contrast, collective fairness expands this notion to include the whole group, emphasizing equitable treatment across all subgroups. It attempts to guarantee that no subgroup suffers disproportionate disadvantages or loss of benefits due to the systemic configuration of the scheme.

Definition 2.2. Collective fairness is realized when the sum of the expected present values (EPV) of all participants’ future benefits is equal to the sum of their initial wealth. Mathematically, this can be represented as

\[ \sum_{i=1}^{L} n_i w_i = \sum_{i=1}^{L}n_iP_0^{i} . \]

Following the methodologies outlined in Chen et al. (Reference Chen, Feng, Wei and Zhao2024), to explicitly quantify actuarial fairness, we introduce the fairness measure $f_{x_i}(s)$ , adapting the concept originally formulated within mutual aid (MA) frameworks into the tontine setting. In the MA context, the fairness measure $f_k$ is defined as the ratio of expected loss to the probability of avoiding that loss event. Analogously, for the tontine framework, $f_{x_i}(s)$ represents the ratio of expected survival benefits to the probability of non-survival. Specifically, we define:

(2.9) \begin{equation} f_{x_i}(s) = \frac{s_{x_i}(s)b_{x_i}(s)}{1 - s_{x_i}(s)}, \end{equation}

where the deterministic survival probability $s_{x_i}(s)$ for cohort $x_i$ at time s is defined by

\[ s_{x_i}(s) = \mathbb{E}[S_{x_i}(s) \mid \mathcal{F}_{t}]= \mathbb{E}\left[\mathbb{E}\left[ \mathbb{I}_{\{\tau_{x_i} \gt s\}} \big| \mathcal{G}_{t+s}^{(i)} \right] \middle| \mathcal{F}_{t}\right]=\mathbb{E}\left[\exp\left\{-\int_0^s \mu_{x_i}(r) dr\right\}\middle| \mathcal{F}_{t}\right]. \]

Here, $i=1,2,\ldots,L$ and $s_{x_i}(s)$ differs from the previously introduced survival probability $S_{x_i}(s)$ through explicit conditioning on the filtration $\mathcal{G}_{t+s}^{(i)}$ , incorporating all historical mortality information except individual-specific death events. The f-value thus indicates the actuarial benefit magnitude for a given cohort i, offering views into the fairness of the benefit distribution.

3.1. Optimal tontine analysis under a homogeneous cohort

In this subsection, we analyze the optimal tontine structure within the context of a homogeneous cohort and Volterra mortality. Considering a homogeneous group composed of n economically similar individuals, each possessing a utility function U, we set the entry age into the tontine as x (with t denote the time of entrance). Operating within the framework of a Volterra mortality model, we employ the individual actuarial fairness principle formalized in Definition 2.1. This principle is combined with the rational behavior of group members aimed at maximizing their expected utility. Guided by this methodology, we derive the OTP structure $ d^*(s) $ for this homogeneous group. Let $\rho$ be the subjective discount factor, r represents the risk free rate, and $w_j$ represents the initial wealth of each individual in the homogeneous cohort and participation rate $\pi_j=1$ . Then our optimization problem is designed as follows:

(3.1) \begin{align} \max_{d(s)} \int_0^\infty e^{-\rho s} \mathbb{E} \left[ U \left( \frac{n w_j d(s)}{N(s)} \right) \mathbb{I}_{\{\tau_{x}\gt s\}} \middle| \mathcal{F}_t \right] ds\notag \\[6pt] \text{subject to} \quad P_0^{j} = \int_0^\infty e^{-r s} \mathbb{E}\left[\frac{n w_j d(s)}{N(s)} \mathbb{I}_{\{\tau_{x}\gt s\}} \middle| \mathcal{F}_t \right]ds = w_j.\end{align}

To formulate and solve the constrained optimization problem, we introduce the Lagrangian, which combines the objective function and the constraints into a single function that can be optimized. The Lagrangian $ \mathcal{L} $ for our problem is defined as

(3.2) \begin{align} \mathcal{L}(d(s), \lambda) = \int_0^\infty e^{-\rho s} \mathbb{E} \left[ S_x(s) \mathbb{E} \left[ U \left( \frac{n w_j d(s)}{N(s)} \right) \middle| \mathcal{G}_{t+s}\right] \middle|\mathcal{F}_{t} \right] \, ds + \notag \\[6pt] \lambda \left( w_j-\int_0^\infty e^{-r s} \mathbb{E} \left[ S_x(s)\mathbb{E} \left[ \frac{n w_j d(s) }{N(s)} \middle| \mathcal{G}_{t+s}\right]\middle|\mathcal{F}_{t}\right] ds \right), \end{align}

where $ \lambda $ is the Lagrange multiplier associated with the wealth constraint.

Next, we need to derive first-order conditions (FOCs):

(3.3) \begin{align}&\mathbb{E} \left[ e^{-\rho s} S_x(s) \mathbb{E} \left[ U' \left( \frac{n w_j d(s)}{N(s)} \right) \frac{n w_j}{N(s)} \middle| \mathcal{G}_{t+s}\right]\middle|\mathcal{F}_{t}\right] \notag \\[6pt] &= \lambda e^{-r s} \mathbb{E} \left[ S_x(s)\mathbb{E} \left[ \frac{n w_j }{N(s)} \middle| \mathcal{G}_{t+s}\right]\middle|\mathcal{F}_{t}\right]. \end{align}

We investigate power and log utility functions within a tontine structure, modeling mortality rates using Volterra processes. Below we provide some simplified representations of the above optimization problem.

Theorem 3.1. Assume a homogeneous group of individuals, each with a power utility function $ U(c) = \frac{c^{1-\gamma}}{1-\gamma} $ , where c denotes consumption and $ \gamma $ is the coefficient of relative risk aversion $(\gamma \neq 1)$ . Given that the mortality rates are influenced by random shocks and modeled by Volterra processes, the OTP structure $ d^*(s) $ for an initial age x and wealth contribution $ w_j $ of each individual is given by

(3.4) \begin{equation} d^*(s) = \frac{1}{w_j } \left( \frac{1}{\lambda}\frac{e^{-(\rho-r)s}\sum_{k=1}^n \binom{n}{k}(\frac{k}{n})^{\gamma}G_{n,k}(s) }{H_{n,x}(s)} \right)^{\frac{1}{\gamma}}, \end{equation}

where $\lambda$ is chosen such that the constraint in (3.1) holds:

\[ \lambda^{\frac{1}{\gamma}} \,:\!=\, \frac{1}{w_j} \int_0^{x^*-x} e^{-(r+\frac{\rho-r}{\gamma})s} H_{n,x}(s)^{-\frac{1}{\gamma}} \left( \sum_{k=1}^n \binom{n}{k} \left(\frac{k}{n}\right)^{\gamma} G_{n,k}(s) \right) ^{\frac{1}{\gamma}} \, ds. \]

\[ H_{n,x}(s) \,:\!=\, \sum_{\ell=1}^n \binom{n}{\ell} ({-}1)^{\ell+1} \exp\left\{-\ell\left(\int_{t}^{t+s} m_x(r)dr-\int_{0}^{t}\eta Y_x(r) dr\right)\right\}\exp\{Z_{x,t}^{(\ell)}(s)\}, \]

\begin{align*} G_{n,k}(s) \,:\!=\, &\sum_{\ell=1}^{n-k} \binom{n-k}{\ell} ({-}1)^\ell\\ & \times \exp\left\{-(\ell+k) \left(\int_{t}^{t+s} m_x(r) \, dr - \int_{0}^{t} \eta Y_x(r) \, dr \right) \right\} \exp\{Z_{x,t}^{(\ell+k)}(s)\}. \end{align*}

For each $ \ell $ in the set $\{1, 2, \ldots, n\}$ , the process $ \{Z_{x,t}^{(\ell)}(s), t \geq 0\} $ is defined with a fixed $ T = x + s $ and assuming $\varphi_\ell (r) \in L^2([0,T])$ , where s ranges from 0 to $ x^* - x $ and $ x^* $ denotes the given terminal age. The process $ Z_{x,t}^{(\ell)}$ are specified by the following stochastic Volterra integral equation:

(3.5) \begin{equation} \left\{\begin{array}{l}\displaystyle Z_{x,t}^{(\ell)} = Z_{x,0}^{(\ell)} + \int_0^{t} \varphi_\ell(T - r) \sigma(Y_{x}^{(\ell)}(r)) \, dW_r - \frac{1}{2} \int_0^t \varphi_\ell(T - r)^2 a(Y^{(\ell)}(r)) \, dr, \\[15pt]\displaystyle Z_{x,0}^{(\ell)} = \int_0^T \left( -\ell \eta + \varphi_\ell(r) b(Y_{x}^{(\ell)}(0)) + \frac{1}{2} \varphi_\ell^2(r) a(Y_{x}^{(\ell)}(0)) \right) dr, \end{array}\right.\end{equation}

where $\varphi_\ell (r)$ solves the Riccati–Volterra equation (RVE)

\[\varphi_\ell (r)= \left( -\ell \eta + b_1 \varphi_\ell(r) + \frac{1}{2} a_1 \varphi_\ell^2(r) \right) K(r), \]

where $\eta$ , $a_1$ , $b_1 \in \mathbb{R}$ are constants and K is a kernel function.

Theorem 3.2. Assume a homogeneous cohort, each with a logarithmic utility function $ U(c) = \ln(c) $ . Considering that mortality rates are subject to random shocks modeled by Volterra processes, the OTP structure $ d^*(s) $ for an initial age x and wealth contribution $ w_j $ of each individual can be calculated as follows:

(3.6) \begin{equation} d^*(s) = \frac{1}{w_j\lambda}\frac{e^{-(\rho-r)s}\mathbb{E}[S_x(s) \mid\mathcal{F}_t]}{H_{n,x}(s)}, \end{equation}

where the Lagrangian multiplier is given by

\[ \lambda = \frac{1}{w_j} \int_0^{x^*-x} e^{-(\rho + r)s/2} \, H_{n,x}(s)^{-1} \, \mathbb{E}[S_x(s) \mid \mathcal{F}_t] \, ds. \]

where $H_{n,x}(s)$ is defined in Theorem 3.1 and

\[ \mathbb{E}[S_x(s) \mid \mathcal{F}_t] = \mathbb{E}[e^{-\int_{t}^{t+s}\mu_x(r) dr}]=\exp\left\{-\left(\int_{t}^{t+s} m_x(r)dr-\int_{0}^{t}\eta Y_x(r) dr\right)\right\}\exp\{Z_{x,t}^{(\ell)}(s)\}, \]

where $\exp\{Z_{x,t}^{(\ell)}(s)\}$ is defined in equation (3.5).

The above theorems provide outcomes for OTPs within a homogeneous cohort. However, in real-world applications, participants often differ in age and wealth, which introduces additional layers of complexity regarding fairness. Next, construct a tontine model that fulfills collective fairness among participants in age-heterogeneous cohorts, addressing the social planner’s problem.

3.2. Multi-cohort tontine design with Volterra mortality

In this subsection, we analyze the optimal design of multi-cohort retirement tontines under the Volterra mortality model, explicitly incorporating the LRD effect of historical mortality deviations. We first establish a general multi-cohort optimization model for heterogeneous groups with different and initial wealth and subsequently study the double-cohort case, where we consider collective fairness introduced in Chen and Rach (Reference Chen and Rach2023).

Consider a tontine pool consisting of L heterogeneous cohorts. Each cohort i (with entry time is t) comprises $n_i$ members who join at age $x_i$ , each contributing an initial wealth $w_i$ . We model the remaining lifetimes of members within each cohort, evaluating their survival probabilities both before and after time s. Our objective is to design a withdrawal rate structure d(s) that not only maximizes the expected utility for each group member but also ensures a fair distribution of wealth between cohorts. We construct the following optimization problem:

(3.7) \begin{align} \max_{d(s)} \int_0^\infty \mathbb{E} \left[\sum_{i=1}^L n_i \beta_i e^{-\rho_i s} U_i \left( \frac{w \pi_i d(s)}{\sum_{l=1}^L N_l(s)} \right) \mathbb{I}_{\{T_{x_i}\gt s\}}\middle| \mathcal{F}_t \right]ds \notag \\[6pt] \text{subject to} \quad \sum_{i=1}^L n_i P_0^i = w = \sum_{i=1}^{L} n_i w_i, \end{align}

where $ U_i $ is the utility functions for all individuals and $\rho_i$ be the subjective discount factor of individual in groups i, assuming homogeneity within each cohort. Furthermore, $\beta_1,\ldots,\beta_L$ be nonnegative numbers adding up to 1 and $ \beta_i = \frac{\bar{w}^{\gamma_{i}}}{\sum_{j=1}^L \bar{w}^{\gamma_{j}}}, \quad \bar{w} = \frac{1}{n} \sum_{i=1}^n w_i. $ To formulate and solve the constrained optimization problem, the Lagrangian $ \mathcal{L} $ for our problem is defined as

(3.8) \begin{align}\mathcal{L}(d(s), \lambda) &= \int_{0}^{\infty} \sum_{i=1}^{L} n_i \beta_i e^{-\rho_i s} \mathbb{E} \left[ S_{x_i}(s) \mathbb{E} \left[ U_i \left( \frac{w \pi_i d(s)}{N(s)} \right) \middle| \mathcal{G}_{t+s}\right]\middle|\mathcal{F}_{t}\right] \, ds \notag \\[6pt]&\quad + \lambda \left( w - \int_{0}^{\infty} \sum_{i=1}^{L} n_i e^{-rs} \mathbb{E} \left[ S_{x_i}(s) \mathbb{E} \left[ \frac{w \pi_i d(s)}{N(s)} \middle| \mathcal{G}_{t+s}\right]\middle|\mathcal{F}_{t}\right] ds\right), \end{align}

where $ \lambda $ is the Lagrange multiplier associated with the wealth constraint.

Next, we need to derive FOCs:

(3.9) \begin{align} &\sum_{i=1}^{L} n_i \beta_i e^{-\rho_i s} \mathbb{E} \left[ S_{x_i}(s) \mathbb{E} \left[ U_i' \left( \frac{w \pi_i d(s)}{\sum_{l=1}^L N_l(s)} \right) \frac{w \pi_i}{\sum_{l=1}^L N_l(s)} \middle| \mathcal{G}_{t+s}\right]\middle|\mathcal{F}_{t}\right] & \notag \\[6pt] &= \lambda e^{-r s} \sum_{i=1}^{L} \left( n_i \pi_i \mathbb{E} \left[S_{x_i}(s) \mathbb{E}\left[\frac{w}{\sum_{l=1}^L N_l(t)} \middle| \mathcal{G}_{t+s}\right]\middle|\mathcal{F}_{t}\right] \right) & \notag \\[6pt] &= \lambda e^{-r s} \sum_{i=1}^{L} \left( n_i \pi_i \mathbb{E} \left[S_{x_i}(s) \sum_{k=0}^{n-1}\frac{w}{k+1} f_{n-1}(k;\, S_{x_1}(s), \ldots, S_{x_L}(s))\middle| \mathcal{F}_{t}\right] \right), \end{align}

where $f_{n-1}(k;\, S_{x_1}(s), \ldots, S_{x_L}(s))$ denotes the PMF of a Poisson binomial distribution for $n-1$ trials with individual survival probabilities $(S_{x_1}(s), \ldots, S_{x_L}(s))$ . This formulation is consistent with the framework established by Fernandez and Williams (2010), which provides a detailed derivation of the distribution.

Formally, consider L heterogeneous cohorts where each cohort i contains $n_i$ initial members with age-specific survival probability $S_{x_i}(s)$ . From the perspective of a specific policyholder in cohort i, the number of surviving members N(s) excluding oneself follows a Poisson binomial distribution with PMF:

\[ f_{n-1}(k;\, S_{x_1}(s), \ldots, S_{x_L}(s))=\Pr\{N(s) = k\} = \sum_{A \in \mathcal{F}_k} \left( \prod_{i \in A} S_{x_i}(s) \right) \left( \prod_{i \notin A} \left(1 - S_{x_i}(s)\right) \right), \]

where $\mathcal{F}_k$ is the collection of all k-element subsets from the $n-1$ members (total pool size $n = \sum_{i=1}^L n_i - 1$ ). Each product term $\prod_{i \in A} S_{x_i}(s)$ quantifies the joint survival probability of exactly k members, while $\prod_{i \notin A} (1 - S_{x_i}(s))$ captures the joint mortality probability of the remaining $n - k - 1$ members.

In the case of power utility, $u(c) = \frac{c^{1-\gamma_i}}{1-\gamma_i},\quad \gamma_i \neq 1$ , one may solve the OTPs of the following form:

\[d^{*}(s)=\left(\frac{\sum_{i=1}^{L}n_i \beta_i e^{-\rho_i s} \mathbb{E}\left[S_{x_i}(s)\sum_{k=0}^{n-1}\left(\frac{w\pi_i}{k+1}\right)^{1-\gamma_i}f_{n-1}(k;\, S_{x_1}(s), \ldots, S_{x_L}(s)) \middle| \mathcal{F}_t\right]}{\lambda e^{-r s} \sum_{i=1}^{L}n_i \pi_i \mathbb{E}\left[S_{x_i}(s)\sum_{k=0}^{n-1}\frac{w}{k+1}f_{n-1}(k;\, S_{x_1}(s), \ldots, S_{x_L}(s)) \middle| \mathcal{F}_t\right]}\right)^{\frac{1}{\gamma_i}}\]

which can be explicitly solved via similar techniques outlined in Theorem 3.1, following from the direct application of results in Abi Jaber et al. (Reference Abi Jaber, Larsson and Pulido2019b).

Although the above formulation provides a theoretical solution for optimal retirement tontines with multiple heterogeneous cohorts under power utility, implementing such a multi-cohort framework in practice poses significant analytical and computational challenges. Specifically, the complexity arises from evaluating high-dimensional conditional expectations involving Poisson binomial distribution $ f_{n-1}(k;\, S_{x_1}(s), \ldots, S_{x_L}(s)) $ . This PMF explicitly depends on each cohort’s survival probabilities, which themselves are random variables due to the historical path-dependent nature of the Volterra mortality model. Consequently, closed-form analytical solutions for optimal payout functions become intractable in most cases, this can be handled in specific forms, via appropriate numerical techniques by considering specific forms of the Volterra mortality model.

Moreover, this model set-up also implicitly allows for different historical mortality experiences of the respective cohorts in the tontine pool. This capability enabled by the Volterra mortality framework distinguishes itself through two mechanisms. First, a deterministic benchmark mortality law establishes baseline survival probabilities that differentiate by age. Second, cohort-specific Volterra processes capture historical deviations from this baseline, thereby reflecting the asymmetrical influences of mortality evolution on optimal payouts. This can be further illustrated from the results under a double-cohort example based on our previous formulation.

Consider again objective and constraint in (3.7) while allowing for tuning parameters $\pi_1,\pi_2$ which differentiates the payouts of two different homogeneous cohorts in a tontine pool, we introduce the Volterral mortality structure as follows. We assume two-dimensional independent Brownian motions ${(W_1(s), W_2(s) : 0 \leq s \leq x^*)}$ driving the deviation of mortality rates from the same benchmark law for the two cohorts. Adopting similar notations as before, from the perspective of any surviving policyholder from cohort 1 at time t after entry:

\begin{align} N_1(s) - 1 \mid \mathcal{G}_{t+s} &\sim \text{Bin}(n_1-1, S_{x_1}(s)), \notag \\[5pt] N_2(s) \mid \mathcal{G}_{t+s} &\sim \text{Bin}(n_2, S_{x_2}(s)).\notag \end{align}

where $N_1(s),N_2(s)$ denote the alive pool members in the cohorts based on Volterra mortality with $n_1, n_2$ being the initial cohort sizes and

\begin{align}S_{x_1}(s) &\,:\!=\, \exp\left\{-\int_{0}^{s} \mu_{x_1} (r) \, dr\right\}, \notag \\[5pt] S_{x_2}(s) &\,:\!=\, \exp\left\{-\int_{0}^{s} \mu_{x_2} (r) \, dr\right\},\notag \end{align}

where the mortality intensities $\{\mu_{x_i}(s) : 0 \leq s \leq x^* - x_i\}$ follow an affine Volterra model ( $i=1,2$ ). Cohorts enter the tontine pool at initial ages $x_1, x_2$ with common terminal age $x^*$ . Similar to the one in the homogeneous cohort case, both of the notional“survival probabilites”above are random variables themselves without conditioning on realized trajectories of the mortality rates up to time $t+s$ . Now it becomes clear that $N(s) \mid \{\mathcal{G}_{t+s}^{(1)} \vee \mathcal{G}_{t+s}^{(2)}\}$ a special case of a general Poisson binomial random variable from the convolution of two distinct binomial distributions with different survival rates. The PM again adopts some general form $ f_{n-1}(k;\, S_{x_1}(s), S_{x_2}(s)) $ outlined in Fernandez and Williams (2010).

In this case of two distinct cohorts, we rewrite the optimization problem in (3.7) as

(3.10) \begin{align} \max_{d(s)} \int_0^\infty \mathbb{E}\left[\left( n_1 e^{-\rho_1 s} \beta_1 U_1 \left(\frac{(n_1 w_1 + n_2 w_2) \pi_1 d(s)}{N_1(s)+N_2(s)}\right) \mathbb{I}_{\{\tau_{x_1} \gt s\}} \right. \right. \nonumber\\[6pt] \left. \left. + n_2 e^{-\rho_2 s} \beta_2 U_2 \left(\frac{(n_1 w_1 + n_2 w_2) \pi_2 d(t)}{N_1(s)+N_2(s)}\right) \mathbb{I}_{\{\tau_{x_2} \gt s\}} \right) \middle| \mathcal{F}_t \right] ds\nonumber\\\text{subject to} \quad \sum_{i=1}^{2} n_i P_0^{i} = w = n_1 w_1 + n_2 w_2, \end{align}

where $ U_1 $ and $ U_2 $ are the utility functions for individuals in cohorts 1 and 2, respectively, assuming homogeneity within each cohort. The subjective discount rates $ \rho_1 $ and $ \rho_2 $ are cohort-specific parameters. From the FOCs in (3.10), we need to calculate $ f_{n-1}(k;\, S_{x_1}(s), S_{x_2}(s)) $ denotes the PMF of a poisson binomial distribution for $ n-1 $ trials with individual survival probabilities $ S_{x_1}(s) $ and $ S_{x_2}(s) $ . To compute the probability mass function (PMF) $f_{n-1}(k;\, S_{x_1}(s),S_{x_2}(s))$ for $ k = 0, 1, 2, \ldots, n - 1$ , we use the convolution method. Each participant’s survival can be represented by a Bernoulli random variable. For cohort 1 (age $ x_1 $ ), the survival probability for each of the $ n_1-1 $ participants (excluding the participant under consideration) is $ S_{x_1}(s) $ . For cohort 2 (age $ x_2 $ ), the survival probability for each participant is $ S_{x_2}(s) $ .

The distribution of N(s) is derived from the convolution of survival probabilities, with probability generating function (PGF):

(3.11) \begin{align} P(z)=(1 - S_{x_1}(s) + zS_{x_1}(s) )^{n_1-1}(1 - S_{x_2}(s) + zS_{x_2}(s) )^{n_2}. \end{align}

This expression represents the probability generating function (PGF) of the Poisson binomial distribution. Expanding P(z) and extracting the coefficients of $ z^k $ gives the probabilities $ f_{n-1}(k;\, S_{x_1}(s),S_{x_2}(s)) $ . Fernandez and Williams (2010) have outlined the general practice in solving for the probability masses $ f_{n-1}(k;\, S_{x_1}(s),S_{x_2}(s)) $ :

(3.12) \begin{align} f_{n-1}(k;\, S_{x_1}(s),S_{x_2}(s)) &= \frac{1}{n} \sum_{m=0}^{n-1} \exp\left( -\nu \frac{2\pi}{n} m k \right) \prod_{l=1}^{2} \left[1 - S_{x_l}(t) + S_{x_l}(t) \exp\left( \nu \frac{2\pi}{n} m \right) \right]^{n_l} \nonumber \\[4pt]&= \frac{1}{n} \sum_{m=0}^{n-1} \exp\left( -\nu \frac{2\pi}{n} m k \right) \left[ 1 - S_{x_1}(s) + S_{x_1}(s) \exp\left( \nu \frac{2\pi}{n} m \right) \right]^{n_1} \nonumber \\[4pt]& \times \left[ 1 - S_{x_2}(s) + S_{x_2}(s) \exp\left( \nu \frac{2\pi}{n} m \right) \right]^{n_2}, \end{align}

where $ \nu = \sqrt{-1} $ is the imaginary unit and $ \frac{2\pi}{n} $ is the angular frequency for discrete Fourier transform (DFT). $ n = n_1 + n_2 $ is the total initial cohort size. This formula leverages Fourier analysis to aggregate individual survival probabilities into the collective distribution.

Theorem 3.3. Consider two heterogeneous cohorts indexed by $i\in\{1,2\}$ , characterized by distinct ages $ x_i $ , subjective discount rates $ \rho_i $ , initial wealth allocation parameters $ \pi_i $ , and cohort-specific relative risk aversion parameters $ \gamma_i \neq 1 $ , sharing a power utility function $ u_i(c) = \frac{c^{1-\gamma_i}}{1-\gamma_i} $ . Assume the survival probability of each cohort follows an affine Volterra mortality model defined as $ S_{x_i}(s) $ . Under these conditions, the OTP function at time s for the two heterogeneous cohorts is given explicitly by

(3.13) \begin{equation} d^*(s) = \left( \frac{1}{\lambda} \frac{\sum_{i=1}^2 \beta_i e^{-(\rho_i - r)s} (\pi_i)^{1-\gamma_i} G_{n_i}(s)}{\sum_{i=1}^2 n_i \pi_i H_{n_i,x_i}(s)} \right)^{\frac{1}{\gamma_i}}. \end{equation}

where

\begin{align} & G_{n_i}(s) \,:\!=\, \mathbb{E}\left[ S_{x_i}(s) \sum_{k=0}^{n-1} \frac{(w\pi_i)^{1-\gamma_i}}{(k+1)^{1-\gamma_i}} f_{n-1}(k;\, S_{x_1}(s),S_{x_2}(s)) \middle| \mathcal{F}_t \right], \notag\\[5pt] & \quad H_{n_i,x_i}(s) \,:\!=\, \mathbb{E}\left[ S_{x_i}(s) \sum_{k=0}^{n-1} \frac{w}{k+1} f_{n-1}(k;\, S_{x_1}(s),S_{x_2}(s)) \middle| \mathcal{F}_t \right].\notag \end{align}

The PMF $f_{n-1}(k;\, S_{x_1}(s),S_{x_2}(s))$ is given as in (3.9), and $ G_{n_i}(s)$ and $H_{n_i,x_i}(s)$ can be separately evaluated given the expected survival probabilities:

(3.14) \begin{align} \mathbb{E}\left[(S_{x_1}(s))^{l} \mid \mathcal{F}_t\right] &= \exp\left\{-l\left(\int_{x_1}^{x_1+s} m_{x_1}(r) \, dr - \int_{0}^{x_1} \eta Y_{x_1}(r) \, dr\right)\right\} \notag\\[5pt]& \times \mathbb{E}\left[\exp\left\{-\ell \int_{0}^{x_1+s} \eta Y_{x_1}(r) \, dr\right\} \, \middle| \, \mathcal{F}_{t}\right], \notag\\[5pt]\mathbb{E}\left[(S_{x_2}(s))^{m} \mid \mathcal{F}_t\right] &= \exp\left\{-m\left(\int_{x_2}^{x_2+s} m_{x_2}(r) \, dr - \int_{0}^{x_2} \eta Y_{x_2}(r) \, dr\right)\right\} \notag\\[5pt]& \times \mathbb{E}\left[\exp\left\{-\ell\int_{0}^{x_2+s} \eta Y_{x_2}(r) \, dr\right\} \, \middle| \, \mathcal{F}_{t}\right].\end{align}

The exponential affine transforms for $i=1,2$ , satisfy $\mathbb{E}[e^{-\int_{0}^{x_i+s}\ell \eta Y_{x_i}(s)ds}|\mathcal{F}_t]=\exp\{Z_{x_i,t}^{\ell}\}$ under which:

(3.15) \begin{equation}\begin{cases} \displaystyle Z_{x_i,t }^{(\ell)} = Z_{x_i, t0}^{(\ell)} + \int_0^{x_i} \varphi_k(T - r) \sigma(Y_{x_i}^{(k)}(r)) \, dW_r - \frac{1}{2} \int_0^t \varphi_\ell(T - r)^2 a(Y_{x_i}(r)) \, dr, \\[2ex] \displaystyle Z_{x_i, 0}^{(\ell)} = \int_0^T \left( -\ell\eta + \varphi_\ell(r) b(Y_{x_i}^{(i)}(0)) + \frac{1}{2} \varphi_\ell^2(r) a(Y_{x_i}^{(i)}(0)) \right) dr, \end{cases} \end{equation}

where $\varphi_\ell(r)$ solves the Riccati–Volterra equation:

\[\varphi_\ell(r) = \left( -\ell\eta + b_1 \varphi_\ell(r) + \frac{1}{2} a_1 \varphi_\ell^2(r) \right) K(r).\]

Theorem 3.4. Assume two age-heterogeneous cohorts of individuals, each with a log utility function $ U(c) = ln(c) $ . Assume that each cohort’s survival probability follows an affine Volterra mortality model $S_{x_i}(t)$ , and that the risk-free rate r is constant. Further assume that the tontine pool is actuarially fair under these mortality and discounting assumptions. Then, the OTP function $d^*(t)$ takes the form

(3.16) \begin{align} d^*(s) = \frac{1}{\lambda} \left( \frac{\sum_{i=1}^2 n_i \beta_i e^{-(\rho_i - r)s} \mathbb{E}\left[ S_{x_i}(s) \sum_{k=0}^{n-1} f_{n-1}(k;\, S_{x_1}(s),S_{x_2}(s)) \mid \mathcal{F}_t ) \right]}{\sum_{i=1}^2 n_i \pi_i H_{n_i,x_i}(s)} \right). \end{align}

where $\lambda$ is the Lagrange multiplier enforcing the actuarial fairness constraint.

In Subsection 2.3, we introduced the fairness measure known as the f-value, which works as an indicator of actuarial fairness among different cohorts in a tontine scheme. In Theorems 3.3 and 3.4, we get the expression for the OTP $d^*(s)$ , then we substituting the expression for $d^*(s)$ into equation (2.9), we get:

\[f_{x_i}(s)=\frac{\pi_i s_{x_i}(s)I_{\{\tau_{x_i}\gt s\}}\sum_{i=1}^{n_i} w_i d^*(s)}{(1-s_{x_i}(s))(N(s)+1)}.\]

To summarize, this section gives the optimal tontine under both power utility and logarithmic utility functions, incorporating the Volterra mortality hypothesis, for homogeneous and age-heterogeneous tontine, respectively. The next section will present numerical simulations and analyses of these theoretical results, validating the model and examining how various parameters influence the optimal strategy. Moreover, we apply the quantitative f -value measure to assess the relative advantages for different groups within heterogeneous cohorts.

4.1. Volterra Vasicek model and simulation process

In the subsection, we introduce the Volterra-type Vasicek mortality model, based on the Gompertz mortality assumption, with the Y process being defined as an affine process. In our numerical analysis of the OTP, the mortality process $ Y_s $ satisfies the following stochastic Volterra integral equation:

(4.1) \begin{equation} Y_s = Y_0 + \int_0^s \frac{(s - r)^{\alpha-1}}{\Gamma(\alpha)} \lambda(\theta - Y_r) dr + \int_0^s \frac{(s - r)^{\alpha-1}}{\Gamma(\alpha)} \sigma dW_r, \quad \end{equation}

where K(t) is the kernel function defined as

(4.2) \begin{equation} K(t) = \frac{t^{\alpha - 1}}{\Gamma(\alpha)}, \quad \alpha \gt 1, \end{equation}

with $\Gamma(\alpha)$ being the Gamma function. This fractional kernel introduces memory effects into the process, allowing for LRD in mortality rates. Research on fractional kernels in rough volatility models has been well documented (Euch and Rosenbaum, 2018; Abi Jaber, et al., Reference Abi Jaber, Bouchard and Illand2019a; Wang et al., Reference Wang, Chiu and Wong2021). This is a fractional VV process with $b_0 =\lambda \theta, b_1 = -\lambda, a_0 = \sigma^2$ . This process belongs to the broad class of Volterra OU processes, and thus key results from Abi Jaber et al. (Reference Abi Jaber, Larsson and Pulido2019b) can be applied. In particular, under the VV model, the function $ \varphi_\ell(s) $ defined in Theorem 3.1 is given as

\[ \varphi_\ell (s) = \int_0^s -\eta \ell E_b(r) dr, \]

where $ E_b(s) $ is specified through the resolvent $ R_b $ of $ \lambda K(t) $ :

\[ E_b(s) = s^{\alpha-1} \sum_{k=0}^{\infty} \frac{({-}\lambda s^{\alpha})^k}{\Gamma(\alpha(k+1))}, \quad R_b(s) = \lambda s^{\alpha-1} \sum_{k=0}^{\infty} \frac{({-}\lambda s^{\alpha})^k}{\Gamma(\alpha(k+1))}. \]

By the conditional expectation of Y, we can explicitly compute $ \exp{\{Z^{(\ell)}_{x,t}(s)\}} $ in $H_{n,x}(s$ ). The specific calculation process is shown in the Appendix A.1. This calculation directly applies the general results from Abi Jaber et al. (Reference Abi Jaber, Larsson and Pulido2019b) and is implemented by recognizing that the resolvent of the fractional kernel K is a Mittag–Leffler function (see Richard et al., Reference Richard, Tan and Yang2021).

For numerical simulations, we discretize the continuous-time stochastic Volterra integral equation using an Euler-type method with fixed step size $ \Delta s = 0.05 $ . Specifically, we define discrete time points $ s_k = k \Delta s $ , for $ k = 0, 1, \ldots, n $ , and approximate the process $ Y_s $ by the discrete approximation $ Y_{s_k} $ . The Euler discretization scheme is given explicitly as follows:

\[ Y_{s_{k+1}} = Y_0 + \sum_{i=0}^{k}\frac{(s_{k+1} - s_i)^{\alpha - 1}}{\Gamma(\alpha)}\,\lambda(\theta - Y_{s_i})\,\Delta s + \sum_{i=0}^{k}\frac{(s_{k+1} - s_i)^{\alpha - 1}}{\Gamma(\alpha)}\,\sigma\,\Delta W, \]

where the increments of the Brownian motion $\Delta W$ are defined as independent over the same grid.

Based on parameters commonly used in the literature, we choose parameters inspired by the work of Wang et al. (Reference Wang, Chiu and Wong2021), and use similar magnitudes. The simulation of the historical trajectory stops at the age of 65 or 70, corresponding to the age when participants enter the tontine scheme.

Table 1. Simulation parameters.

Figure 1. 3 sample paths of the specified VV process under 65.

Figure 1 illustrates three sample trajectories generated under the Volterra mortality model to simulate the historical mortality experience of a cohort aged 65. These trajectories play a significant role in calculating the optimal tontine. We emphasize the significant effect of the long-memory properties inherent in affine Volterra processes on the structure of fair tontine payouts and demonstrate that comparing these results with the benchmark Gompertz mortality model allows for a more precise analysis of payment differences across different cohorts. While the existing literature offers a relatively well-recognized framework for tontine fairness, there has yet to be an in-depth exploration of the long range dependence in mortality structures.

This simulation process establishes a foundation for computing the OTPs based on the generated mortality trajectories, which we discuss in the next subsection. In the subsequent analysis, we will select one sample path as the basis for calculating the optimal tontine.

4.2. Optimal tontine payout

In this subsection, we present the numerical results for the optimal tontine payout (OTP) based on the VV mortality model, analyzing its performance across different population groups and utility functions. To compute the OTP, we use the mortality trajectories generated in Subsection 4.1 and numerically solve the continuous-time stochastic differential equation using the discrete-time Euler method. Similar to Chen and Rach (Reference Chen and Rach2023), we apply the Gompertz law for the deterministic counterpart m(s) for the mortality rate $\mu$ . This is required as we specify $\eta Y$ to be the random deviation from m that fluctuates with time. We have chosen the set of parameters based on the empirical study conducted by Milevsky (Reference Milevsky2020) on calibrating the Gompertz law to different populations.

Table 2. Base case parameters setup. We assume individuals belong to cohorts $i=1,2$ .

Figure 2. Illustration of the Volterra–Vasicek mortality process and the resulting survival probabilities.

We first select a representative mortality trajectory and calculate the OTP along this path, as shown in Figure 2(a). The results from the conditional expectation of the historical mortality trajectory, demonstrate significant path dependence, reflecting the long-memory property of the VV model. This single-path realization demonstrates significant fluctuations around zero, consistent with the long-memory property intrinsic to the fractional Volterra specification.

However, single-path trajectories alone cannot fully elucidate LRD’s cumulative impact on survival outcomes. To address this, Figure 2(b) presents a comparison of survival probabilities under three distinct mortality assumptions: the traditional Gompertz model, the affine Vasicek model without LRD, and the Volterra–Vasicek model incorporating LRD. While the Gompertz and affine Vasicek curves show close alignment, highlighting their functional similarity in the absence of long-memory effects. The Volterra–Vasicek survival curve initially falls below both the Gompertz and affine curves, reflecting a faster early decline; however, as historical mortality dependence accumulates, it decays more slowly and ultimately produces higher long-term survival probabilities.

For a homogeneous population, Figure 3 illustrates the comparative dynamics of the OTP over time under three distinct mortality frameworks: the traditional Gompertz model (with random shock $\varepsilon$ ), the affine mortality model and the Volterra mortality model incorporating long-range dependence (LRD). The left panel corresponds to the logarithmic utility function, while the right panel corresponds to the power utility function. These comparisons offer insights into how different mortality assumptions can significantly affect the optimal payout strategies for cohorts with homogeneous age characteristics, highlighting the role of historical mortality dependence in payout design.

Figure 3. Optimal tontine results under Gompertz, Affine and Volterra mortality models in a homogeneous cohort.

Specially, as shown in Figure 3, the optimal tontine results under Gompertz and affine mortality models in a homogeneous cohort, do not exhibit a significant difference. In contrast, the Volterra-based OTP exhibits a sharper initial decline relative to Gompertz. Arguably this is due to the analytical solution is of a form that is different from the ones presented in Chen et al. (Reference Chen, Hieber and Klein2019, Reference Chen, Hieber and Rach2021) due to the inclusion of Volterra processes. Intuitively, this is also partially explained by the LRD structure. As shown in Figure 2(a), the green solid line representing Y lies below the black dashed zero line up to age 65. Deviations of the Y trajectory from zero, with over 50% of the values being negative, indicate that the cohort’s mortality rate before age 65 is lower on average than the mortality rate projected by the reference mortality law.

Such trajectories imply lower-than-expected mortality experiences on average for the cohort before age 65, consequently translating into higher anticipated longevity risks. Due to the model’s dependency on historical mortality paths, the implied Volterra survival probabilities become higher, reflecting increased longevity uncertainty. As a result, an early-stage payout reduction emerges as a necessary risk compensation mechanism, mitigating potential longevity risk through conservative initial payments.

For mixed-age populations, we analyze two distinct age cohorts of tontine participants with collective fairness. Figure 4 extends this analysis to populations characterized by age heterogeneity, explicitly differentiating two distinct age cohorts. Similar to the homogeneous scenario, the left panel again presents results for the logarithmic utility function, whereas the right panel illustrates those under power utility. Each panel presents three mortality specifications, traditional Gompertz, affine Vasicek, and Volterra–Vasicek. This extension facilitates a direct assessment of mortality model choice under collective fairness, affecting payout outcomes for distinct age cohorts.

Figure 4. Optimal tontine results under Gompertz, Affine and Volterra mortality models in two-age heterogenous cohorts.

From Figure 4, the Gompertz and affine curves almost coincide over the entire horizon. While introducing the Volterra kernel generates a markedly steeper initial decline, followed by a smoother taper, in both utility cases. These curves demonstrate the payment decay over time for each strategy, highlighting variations in time sensitivity and risk preferences. The general trend aligns with findings in the existing literature, such as Chen and Rach (Reference Chen and Rach2023).

When comparing the two utilities, we observe that the introduction of affine mortality results in a significant reduction in the initial optimal tontine, with a notably steeper slope, followed by a gradual and smoother decrease over time. This suggests that the incorporation of the Y process, accounting for historical factors such as past mortality levels leads the model to predict greater future longevity risk, prompting a more conservative payout in the early phases. This model structure more accurately reflects the real-life uncertainty surrounding longevity risk. It is evident that models incorporating the Volterra process exhibit a faster initial decline in pension payments, likely due to their reliance on historical data, which makes the model’s estimates of longevity risk more conservative. This early decline serves as an insurance against uncertain risks associated with future life extension.

4.3. Fairness analysis

In this subsection, we examine the fairness across different age groups within the tontine structure, using the fairness measure f-value proposed by Feng et al. (Reference Feng, Liu and Zhang2024). The f-value quantifies the actual payments received by each group, accounting for differences in survival probability and mortality rates, thereby revealing fairness tendencies in the payment process. We analyzed tontine participants from younger and older groups, calculating the f-value for each group. The results indicate that younger participants typically receive higher relative gains, reflecting a certain fairness bias. This is mainly due to the higher survival probabilities in the younger group, which allow them to claim a larger share of the tontine pool. Consequently, older participants may face a relative disadvantage under the same payment structure, displaying a degree of unfairness.

Figure 5 plots the actuarial fairness measure (f-value) dynamics for two entry ages, 65 and 70, under three mortality specifications – Gompertz, affine Vasicek (no LRD), and Volterra–Vasicek (LRD). For both age settings (65 and 70), the Gompertz and affine curves are virtually indistinguishable, whereas the Volterra curve departs markedly owing to the presence of long-range dependence. Although actuarial fairness is comprehensively maintained across groups, there remains an inherent fairness issue disproportionately affecting older participants. Both age cohorts exhibit significant initial declines in their respective f -values, but notably, the group entering at age 70 consistently demonstrates lower f -values over most of the evaluation period. This observation indicates that older participants face greater actuarial disadvantages compared to their younger counterparts, highlighting an age-based inequity within the tontine arrangement. This finding aligns with Chen and Rach (Reference Chen and Rach2023), which shows that actuarial fairness optimized for a common cohort payout $ d^*(s) $ may not be optimal for specific subgroups.

Figure 5. Actuarial fairness measures for two heterogeneous age cohorts under distinct mortality models.

The incorporation of the Volterra mortality model, which introduces LRD into mortality rates, influences the f-values by affecting the differences between age groups and intensifies these disparities. The long-memory effects captured by the Volterra process mean that historical mortality experiences have a more pronounced influence on current survival probabilities. In the early to mid term, Volterra–Vasicek yields slightly more conservative survival estimates for older cohorts, but as time advances its decay rate slows and its survival probabilities eventually exceed those of both Gompertz and affine Vasicek models. This behavior produces notably higher initial f-values and a wider inter-cohort spread, while all models converge toward zero over the horizon. Under the Gompertz model initial f-values fall within a moderate range. And the affine Vasicek extension closely mirrors this profile, differing only by small cohort-specific shifts. By contrast, the Volterra–Vasicek curve displays a considerably broader dispersion of initial f-values, extending to values roughly an order of magnitude higher. The application of the f-value fairness measure, adapted from Chen et al. (Reference Chen, Feng, Wei and Zhao2024) for tontines by redefining survival as the triggering event, effectively captures these nuances. Our findings confirm the practical effectiveness of this fairness metric, underscoring the necessity of explicitly considering age heterogeneity and historical mortality dependence in ensuring equitable actuarial outcomes across diverse age cohorts.