2.1. Basic notations

Consider a portfolio of $N_0$ policyholders aged x at time 0. The random variable $I_i(0,k)$ is equal to 1 if policyholder i survives to time $k\geq 0$ , and to 0 otherwise, for $i=1,...,N_0$ . The number of survivors at time k is $N_k=\sum_{i=1}^{N_0}I_i(0,k)$ . For the $N_j$ surviving policyholders at time j, the random variable $I_{l}(\,j,k)$ is equal to 1 if policyholder l survives from time j to time k, and to 0 otherwise, with $j\leq k$ and $l=1,...,N_j$ . Note that the order of policyholders is not relevant in the present context. Note also that the number of policyholders at time k can also be written as $N_k=\sum_{i=1}^{N_j}I_i(\,j,k)$ . Let I(j,k) be the random variable representing the proportion of survivors from time j to time k. Specifically, for $N_j\gt0$ and $j\leq k$ :

\begin{align*}I(\,j,k)=\frac{1}{N_j}\underset{i=1}{\overset{N_j}{\sum}}I_i(\,j,k)=\frac{I(0,k)}{I(0,j)}.\end{align*}

This survival index is only defined for $N_j\gt0$ . If $N_j=0$ , meaning that all policyholders have died by time j, the index I(j,k) becomes irrelevant as no further benefits will be paid.

The remaining lifetimes of policyholders are assumed to be identically distributed and conditionally independent. Conditional independence means that their survival is influenced by shared external factors such as medical advancements or socioeconomic dynamics. The generic notation $\Theta$ , whose distribution is not crucial in the analytical derivations, is used to represent those common external factors. This notation highlights the dependence between the survival indicators and the randomness of the survival probabilities. Namely, the influence of external factors $\Theta$ on the survival probability is highlighted in the notations such that ${}_{k-j}p_{x+j}(\Theta)$ is the $(k-j)$ -survival probability at time j for a given realization of $\Theta$ :

\begin{align*}{}_{k-j}p_{x+j}(\Theta) = \mathbb{P}\left[I_i(\,j,k)=1|\Theta\right] = \mathbb{E}\left[I(\,j,k)|\Theta\right]=1 - {}_{k-j}q_{x+j}(\Theta).\end{align*}

The corresponding ‘unconditional’ survival probabilities are denoted by the conventional actuarial symbol ${}_{k-j}p_{x+j}=\mathbb{E}\left[{}_{k-j}p_{x+j}(\Theta)\right].$ The influence of external factors $\Theta$ on the dependence between policyholders’ remaining lifetimes is such that for $i\neq l$ :

\begin{align*}\mathbb{P}\left[I_i(0,k)= I_l(0,k)=1\right] =\mathbb{E}\left[\mathbb{P}\left[I_i(0,k)=1|\Theta\right]\mathbb{P}\left[I_l(0,k)=1|\Theta\right]\right]=\mathbb{E}\left[{}_kp_x(\Theta)^2\right].\end{align*}

For $j\leq k$ , the random survival probabilities ${}_{k-j}p_{x+j}(\Theta)$ are assumed to be almost surely strictly positive over the time interval [0, T] of interest. Note that for $k\geq j$ , conditionally on $\Theta$ , the property ${}_kp_x(\Theta)={}_jp_x(\Theta) \times {}_{k-j}p_{x+j}(\Theta)$ holds. However, taking the expectation on both sides, ${}_kp_x$ is not necessarily equal to the product of expectations ${}_jp_x \times {}_{k-j}p_{x+j}$ .

Throughout the paper, the term ‘diversifiable risk’ refers to the risk associated with portfolio size at time 0, with larger $N_0$ associated with lower diversifiable risk. The term ‘systematic risk’ refers to the risk stemming from the uncertainty in the survival probabilities, with higher uncertainty (or higher variance in this paper) associated with higher systematic risk.

2.2. Annuities under risk sharing

At time 0, all $N_0$ policyholders underwrite a T-year term annuity that pays an initially agreed-upon benefit amount b in arrears upon survival. At contract inception, each policyholder pays a premium $P_0=(1+\varphi)\pi$ , where $\pi$ is the actuarially fair premium, given by the expected present value of all benefits $\pi=ba_{x:\overline{T}\kern-1.1pt\raise-1.55pt\hbox{$|$}}$ , with $a_{x:\overline{T}\kern-1.1pt\raise-1.55pt\hbox{$|$}}$ denoting the annuity factor, and $\varphi\gt0$ representing the risk loading charged by the insurer.

The annuity is assumed to be priced under the DEP introduced in Hanbali et al. (Reference Hanbali, Denuit, Dhaene and Trufin2019). Under a classical annuity, the yearly shortfall arising from the deviation between the available and the required reserves is entirely borne by the insurer. In contrast, under the DEP, the yearly shortfall is shared with policyholders. This means that at time k, the initially agreed-upon benefit b payable to the $N_k$ surviving policyholders is adjusted by $P_k$ , where the adjustment reflects the deviation caused over the year $(k-1,k)$ . In particular, the amount effectively payable to the $N_k$ surviving policyholders at time k is $b-P_k$ .

For $N_{k}\gt0$ , from Proposition 1 in Hanbali (Reference Hanbali2025), the adjustment $P_k$ is given by

(2.1) \begin{equation}P_k=\alpha_k \left(1+\varphi\right)b\frac{\lambda(0,k)}{v^k}\Delta_k, \quad k=1,...,T,\end{equation}

such that

(2.2) \begin{equation} \Delta_k = \frac{1}{\tilde{p}_{x+k-1}} - \frac{1}{I(k-1,k)}, \end{equation}

and

(2.3) \begin{equation} \lambda(\,j,k) = v^{k+j-1}a_{x+k-1:\overline{T-k+1}\kern-1.1pt\raise-1.55pt\hbox{$|$}} \quad \text{and} \quad \tilde{p}_{x+k-1}=\frac{a_{x+k-1:\overline{T-k+1}\kern-1.1pt\raise-1.55pt\hbox{$|$}}}{v\left(1+a_{x+k:\overline{T-k}\kern-1.1pt\raise-1.55pt\hbox{$|$}}\right)}. \end{equation}

Here, $v=(1+r)^{-1}$ is the yearly discount factor (assumed constant), and $\alpha_k\in[0,1]$ is the share of the shortfall transferred back to policyholders at time k. Note that $I(k-1,k)$ is not defined for $N_{k-1}=0$ , as all policyholders have died by time $k-1$ .

An annuity without adjustments (henceforth, a classical annuity) and risk-loading $\psi$ is obtained by setting $\alpha_1=...=\alpha_T=0$ and $\varphi=\psi$ .

For $\underline{\alpha}=\left(\alpha_1,...,\alpha_T\right)$ , the loss random present value of policyholder i at time 0 is

(2.4) \begin{equation}L_i^{pol}(\underline{\alpha},\varphi)=P_0 - b\underset{k=1}{\overset{T}{\sum}}v^kI_i(0,k)+ (1+\varphi)b\underset{k=1}{\overset{T}{\sum}}\lambda(0,k)I_i(0,k)\alpha_k\Delta_k ,\end{equation}

whereas the per-policy loss random present value of the insurer at time 0 is

(2.5) \begin{equation}L^{ins}(\underline{\alpha},\varphi)= b\underset{k=1}{\overset{T}{\sum}}v^kI(0,k) - P_0 + (1+\varphi)b\frac{1}{N_0}\underset{i=1}{\overset{N_0}{\sum}}\underset{k=1}{\overset{T}{\sum}}\lambda(0,k)I_i(0,k)(1-\alpha_k)\Delta_k.\end{equation}

Specifically, for policyholder i, a premium $P_0$ is paid at time 0, and at each time k, they receive the adjusted benefit $b-P_k$ . From the insurer’s perspective, each of the $N_0$ policyholders pays the premium $P_0$ , and the $N_k$ survivors at time k receive the adjusted benefit $b-P_k$ . On top of that, at time k, the insurer covers the shortfall, which is given by $(1+\varphi)b\frac{\lambda(0,k)}{v^k}\Delta_k$ .

2.3. Problem formulation

This paper studies optimal risk-sharing designs for annuities under risk sharing from the perspectives of both policyholders and the insurer. Specifically, the goal is to determine the initial loading under risk-sharing $\varphi$ , and the yearly shares of transferred risk $\alpha_1,...,\alpha_T$ , which are optimal for both parties. This consists of a two-objective maximization problem, where the objective is to find PO $\underline{\alpha}$ and $\varphi$ within some given joint admissibility set $\mathcal{A}$ .

Both parties are endowed with mean-variance preference functions denoted $\rho^{pol}$ and $\rho^{ins}$ for the policyholder and the insurer. Despite its limitation of focusing only on the first two moments and potentially ignoring tail information, the mean-variance preference function has appealing properties. It quantifies preferences by combining the expected value of the net profit with sensitivity to variability. It also approximates the CRRA utility function. In addition, it allows for closed-form expressions that ease the simulation process, offering insights into how the time-varying risk-sharing can be implemented (see Pulley, Reference Pulley1981; Fahrenwaldt and Sun, Reference Fahrenwaldt and Sun2020).

The objective functions of interest are

\begin{eqnarray*} \rho^{pol}(\underline{\alpha},\varphi) &=& \mathbb{E}\left[w_i -L_i^{pol}(\underline{\alpha},\varphi)\right] - \gamma \mathbb{V}\left[w_i -L_i^{pol}(\underline{\alpha},\varphi)\right],\\ \rho^{ins}(\underline{\alpha},\varphi)&=& \mathbb{E}\left[w -L^{ins}(\underline{\alpha},\varphi)\right] - \delta \mathbb{V}\left[w -L^{ins}(\underline{\alpha},\varphi)\right], \end{eqnarray*}

where $\gamma\gt0$ and $\delta\gt0$ measure the sensitivity of each party to variability, and $w_i$ and w are the initial wealth of policyholder i and the insurer per-policy, respectively.

Unlike in several studies on risk sharing in insurance, for a given set of admissible premium loading and adjusted benefits (or premium and indemnity in non-life insurance), the set of objective functions to be maximized is not concave. As a result, yielding all PO $\underline{\alpha}$ and $\varphi$ cannot be achieved simply by maximizing a linear combination of the objective functions, as in for example Asimit et al. (Reference Asimit, Cheung, Chong and Hu2020). However, solving the optimization problem involving the linear combination of the objective functions will yield a set of PO solutions (Boyd and Vendenberghe, Reference Boyd and Vendenberghe2004). The optimization of the linear combination of $\rho^{pol}(\underline{\alpha},\varphi)$ and $\rho^{ins}(\underline{\alpha},\varphi)$ is therefore the primary focus of this paper. In particular, for each $\nu \in (0,1)$ , the goal is to solve the following optimization problem:

(2.6) \begin{eqnarray}&&\arg \underset{\underline{\alpha},\varphi}{\max} \ \ \mathcal{F}_{\nu}(\underline{\alpha},\varphi),\end{eqnarray}

(2.7) \begin{eqnarray}&&\text{subject to } (\underline{\alpha},\varphi) \in \mathcal{A}. \end{eqnarray}

where the objective function $\mathcal{F}_{\nu}(\underline{\alpha},\varphi)$ to be maximized is given by

\begin{align*}\mathcal{F}_{\nu}(\underline{\alpha},\varphi)=\nu \rho^{pol}(\underline{\alpha},\varphi) + (1-\nu) \rho^{ins}(\underline{\alpha},\varphi).\end{align*}

Regarding the admissible set, first, $\alpha_k\gt 1$ implies that the policyholder covers more than the required shortfall, while $\alpha_k\lt0$ implies that the policyholder is compensated when the insurer experiences a loss. Thus, the natural admissible set for $\underline{\alpha}$ is $[0,1]^T$ . For the risk loading, by definition, $\varphi\geq 0$ . Regarding the upper bound, while one option is to allow the loading to be unbounded, a more consistent approach with the goal of enhancing annuity affordability is to cap $\varphi$ at $\psi$ , where $\psi$ is the loading of a benchmark classical annuity. That is, in counterpart for sharing the risk, policyholders’ time 0 premium is cheaper than that of the classical annuity. Thus, the joint admissible set $\mathcal{A}$ is $[0,1]^T \times [0,\psi]$ , where $\psi$ is a given loading for a classical annuity.

3.1. Alternative expression of the objective function

The following proposition provides expressions for the mean-variance preference functions $\rho^{pol}(\underline{\alpha},\varphi)$ and $ \rho^{ins}(\underline{\alpha},\varphi)$ in function of the distributions of the survival probabilities and portfolio size. These expressions will be useful in the analytical derivations and the numerical study, as they avoid the need for nested simulations of the survival indices and variance reduction techniques. The detailed derivations as well as the expressions for the quantities $\mu$ ’s and $\sigma$ ’s are provided in Appendix A.1.

Proposition 1. The mean-variance preference functions of the policyholder and the insurer under risk sharing are given by

(3.1) \begin{eqnarray}\rho^{pol}(\underline{\alpha},\varphi)&=&w_i - \gamma b^2\sigma_a - b a_{x:\overline{T}\kern-1.1pt\raise-1.55pt\hbox{$|$}} \varphi + b(1+\varphi)\underset{k=1}{\overset{T}{\sum}}\alpha_k\left(2\gamma b\sigma_{p,a}^{(k)}-\mu^{(k,k)} \right) \nonumber\\ &&- \gamma b^2(1+\varphi)^2 \left(\underset{k=1}{\overset{T}{\sum}}\alpha_k^2\sigma_p^{(k,k)} + 2\underset{k=2}{\overset{T}{\sum}}\underset{j=1}{\overset{k-1}{\sum}}\alpha_k\alpha_j\sigma_p^{(\,j,k)}\right),\end{eqnarray}

and

(3.2) \begin{eqnarray} \rho^{ins}(\underline{\alpha},\varphi)&=&w -\delta b^2\tilde{\sigma}_a+ ba_{x:\overline{T}\kern-1.1pt\raise-1.55pt\hbox{$|$}}\varphi - b (1+\varphi)\underset{k=1}{\overset{T}{\sum}}(1-\alpha_k)\left(2\delta b\tilde{\sigma}_{p,a}^{(k)}+\mu^{(k,k)}\right) \nonumber\\ &&-\delta b^2(1+\varphi)^2\left(\underset{k=1}{\overset{T}{\sum}}(1-\alpha_k)^2\tilde{\sigma}_p^{(k,k)} + 2\underset{k=2}{\overset{T}{\sum}}\underset{j=1}{\overset{k-1}{\sum}}(1-\alpha_k)(1-\alpha_j)\tilde{\sigma}_p^{(\,j,k)}\right), \end{eqnarray}

where the expressions of $\mu^{(k,k)}$ , $\sigma_a$ , $\tilde{\sigma}_a$ , $\sigma_p^{(\,j,k)}$ , $\tilde{\sigma}_p^{(\,j,k)}$ , $\sigma_{p,a}^{(k)}$ , and $\tilde{\sigma}_{p,a}^{(k)}$ are given in Appendix A.1.

Proposition 1 shows that the shape of the objective function $\mathcal{F}_{\nu}(\underline{\alpha}, \varphi)$ presents a prohibitive challenge. Specifically, the function is quartic when the parameters $\underline{\alpha}$ and $\varphi$ are optimized jointly. Furthermore, based on the derivations in Appendix A.8, the Hessian matrix is given by

(3.3) \begin{equation}H_{\left(\varphi,\underline{\alpha}\right)}=\left(\begin{array}{c@{\quad}c}\frac{\partial^2 \mathcal{F}_{\nu}(\underline{\alpha},\varphi)}{\partial \varphi^2} & Z_{\left(\varphi,\underline{\alpha}\right)}^{\top}\\Z_{\left(\varphi,\underline{\alpha}\right)} & \Sigma_{\left(\varphi,\underline{\alpha}\right)}\end{array}\right),\end{equation}

where $Z_{\left(\varphi,\underline{\alpha}\right)}\in\mathbb{R}^{T}$ is a vector whose k-th component is $\frac{\partial^2 \mathcal{F}_{\nu}(\underline{\alpha},\varphi)}{\partial \varphi\partial \alpha_k}$ , and $\Sigma_{\left(\varphi,\underline{\alpha}\right)}$ is a TxT symmetric matrix given by $\left(\frac{\partial^2 \mathcal{F}_{\nu}(\underline{\alpha},\varphi)}{\partial \alpha_k\partial \alpha_j}\right)_{j,k}$ . This Hessian matrix depends on the arguments of the problem, where its entries are quadratic functions, making it impractical to analytically establish the range of values leading to negative semi-definiteness (NSD) as a condition for global concavity. This indicates that solving for globally optimal solutions in the present paper is challenging and computationally intensive using conventional optimization techniques. This is consistent with Ahmadi et al. (Reference Ahmadi, Olshevsky, Parrilo and Tsitsiklis2011) who note that establishing the concavity of the Hessian in quartic optimization is itself an NP-hard problem.

The second challenge, albeit less prohibitive numerically, lies in the numerous constraints imposed on the optimization problem. The admissible set $\mathcal{A}$ constrains $\underline{\alpha}$ and $\varphi$ within the bounds of $[0,1]^T \times [0, \psi]$ . Solving the Lagrangian function requires considering all possible combinations of binding and nonbinding constraints for the Lagrange multipliers. For each $\varphi$ and $\alpha_k$ , these multipliers can either be equal to 0 or strictly positive, depending on whether the constraints are active. This results in too many combinations to check, which is manageable for small T, but quickly becomes impractical for larger values of T which are to be expected in realistic annuity contracts.

The present paper addresses the complexity of the Pareto optimization problem in two stages. The first stage involves an intensive numerical approach, where several algorithms are applied to attempt to identify the PO values of $\varphi$ and $\underline{\alpha}$ which solve the optimization problem (2.6)–(2.7). This numerical approach, along with the corresponding results, is outlined in the following subsetion. The second stage leverages patterns observed in the numerical solutions to derive a heuristic analytical solution. The closed-form analytical expression resulting from this heuristic, as well as its efficiency, is examined in the final part of this section.

3.2. Numerical solution

3.2.1. Input

Mortality is modeled using the Lee–Carter model calibrated with Australian data from 1950 to 2018 and for ages 25 to 99 (Lee and Carter, Reference Lee and Carter1992), where the time-dependent factor representing improvements in mortality is modeled as a random walk with drift. Multiple paths of the time components are generated, leading the paths of the central death rates, and eventually to the paths of the yearly survival probabilities $p_{x+k}(\Theta)$ . The latter are multiplied to obtain the k-year survival probabilities ${}_kp_x(\Theta)$ .

To isolate the effect of systematic risk, the baseline simulated yearly probabilities are adjusted by a variance factor, such that:

\begin{align*}p_{x+k}(\Theta^{(s)}) = \left(p_{x+k}(\Theta)- \mathbb{E}\left[p_{x+k}(\Theta)\right]\right)\times \frac{s}{2} + \mathbb{E}\left[p_{x+k}(\Theta)\right].\end{align*}

In particular, the variance of $p_{x+k}(\Theta^{(s)})$ is $\frac{s}{2}$ times larger than the variance of the baseline probabilities. The case study presented here considers $\frac{s}{2}=1$ (i.e., the baseline probabilities), as well as $\frac{s}{2}=3$ (i.e., high systematic risk), for which all probabilities are in the unit interval despite the linear transformation with the variance factor.

The effect of diversifiable risk is controlled by portfolio size, with two sizes considered, $N_0=100$ and $N_0=\infty$ . The case study evaluates the problem with terms $T=5$ , 10, 15, 20, and 25 and two ages $x=55$ and $x=70$ . The risk aversion parameters $\gamma$ and $\delta$ are equal to either $0.5$ , interpreted as low sensitivity to variance, or equal to 5, interpreted as high sensitivity to variance. Interest rate is either equal to 0% or to 5% per annum. The loading of the classical annuity $\psi$ is equal to either $0.1$ or $0.3$ . Throughout, it is assumed that the policy pays $b=1$ per annum.

Only selected results are presented in this numerical study, as most combinations of input parameters lead to similar patterns.

3.2.2. Optimization approach

The combination of a quartic objective function and multiple constraints renders the problem highly challenging to solve. As a result, it is important to ensure that the numerical solution is acceptable, even though it is difficult to guarantee that it is the true global optimum. The solution process in this paper follows two stages: an extensive numerical approach, followed by identifying patterns that lead to a simplified analytical heuristic.

The numerical approach involves applying three optimization methods and selecting the one that provides the highest objective function value as the optimal solution. All three methods attempt to avoid trapping in local optima.

The first method uses diffused starting points within the bounds of $\underline{\alpha}$ and $\varphi$ for the joint optimization, which is solved using Byrd et al. (Reference Byrd, Lu, Nocedal and Zhu1995)’s algorithm. The initial values for this approach include combinations of $\frac{\varphi}{\psi}\in\{ 0, 0.25, 0.5,0.75,1\}$ , paired with different $\alpha_k$ configurations. These configurations are constant risk sharing ( $\alpha_k = \alpha \in \{0, 0.25, 0.5, 0.75, 1\}$ ), increasing risk sharing ( $\alpha_k = \frac{k}{T}$ ), decreasing risk sharing ( $\alpha_k = 1 - \frac{k}{T}$ ), as well as U- and inverted U-shapes of $\alpha_k$ , that is, $\alpha_k$ equals 0 during the first and last quarters of the policy and 1 otherwise, or the reverse.

The second involves a separate optimization of $\varphi$ and $\underline{\alpha}$ , rather than the original joint optimization. Specifically, $\underline{\alpha}$ are optimized for different fixed values of $\varphi$ , after which the best $\varphi$ is selected based on the highest objective function value. The resulting $\varphi$ and $\underline{\alpha}$ are then used as initial values for the joint optimization, again solved using Byrd et al. (Reference Byrd, Lu, Nocedal and Zhu1995)’s algorithm.

The third method evaluates the objective function at several uniformly generated points $(\underline{\alpha},\varphi)\in\mathcal{A}$ and approximates the objective function using OLS as

\begin{align*}\mathcal{F}_\nu(\underline{\alpha},\varphi)\approx \hat{\mathcal{F}}_\nu(\underline{\alpha},\varphi)=\beta_0 + \beta_{\varphi}\varphi + \gamma_{\varphi}\varphi^2 + \sum_{k=1}^T\left(\beta_k\alpha_k+\gamma_k\alpha_k^2\right).\end{align*}

Being quadratic, this approximation has optimal $(\underline{\alpha},\varphi)$ which are substantially easier to identify. Once calculated, they are used as initial points for Byrd et al. (Reference Byrd, Lu, Nocedal and Zhu1995)’s algorithm to optimize the original objective function $\mathcal{F}_\nu$ .

For each combination of input values, the point $\left(\underline{\alpha}, \varphi\right)$ that yields the highest value of the objective function is retained. The next step is to verify whether the Hessian matrix $H_{\left(\underline{\alpha}, \varphi\right)}$ defined in (3.3) is NSD at this point. However, this verification proved challenging. In particular, solving for the eigenvalues $\underline{\lambda}$ resulted in various combinations of signs, with all solutions leading to small values of $|H_{\left(\underline{\alpha}, \varphi\right)}-\lambda I |$ . This issue likely stems from the fact that the determinant of the Hessian matrix itself is often very small.

Nevertheless, for all combinations of input values, it was possible to find all negative eigenvalues of the Hessian matrix evaluated at the optimal $\left(\underline{\alpha}, \varphi\right)$ . Across all input values, the maximum absolute value of the determinant of $|H_{\left(\underline{\alpha}, \varphi\right)}-\lambda I |$ with all negative eigenvalues was approximately $3.8 \times 10^{-10}$ . This result is interpreted as indicative of the Hessian’s NSD at the estimated optimal $\left(\underline{\alpha}, \varphi\right)$ .

3.2.3. Results

Before presenting the main results, an analysis of individually optimal contracts was conducted. Although the results are not reported here, the key findings are as follows: the objectives of the policyholder and insurer are generally in conflict regarding the initial loading $\varphi$ across nearly all input parameters. There is also conflict over the risk-sharing parameter $\alpha_k$ , though in cases with larger portfolios ( $N_0 = \infty$ ), both parties align on the optimal $\alpha_k=1$ for much of the contract duration. Additionally, neither the interest rate nor the initial loading $\psi$ plays a significant role in shaping the optimal designs.

Figure 1 presents the PO risk-sharing designs for policyholders aged 70 years with a policy term of 25 years, under a 0% interest rate and $\psi = 0.1$ . The results are consistent across other interest rates and values of $\psi$ , as well as for age 55. Each row of the figure displays the results for different combinations of the risk aversion parameters $\gamma$ and $\delta$ . The left and middle columns show the PO $\alpha_k$ for portfolio sizes $N_0 = 100$ and $N_0 = \infty$ , where lighter shades correspond to more weight on the policyholder’s objective ( $\nu$ close to 1), and darker shades reflect more weight on the insurer’s objective ( $\nu$ close to 0). The red curves correspond to the baseline systematic risk case, and the blue ones correspond to the high systematic risk case. The right column displays the values of $\frac{\varphi}{\psi}$ as a function of $\nu$ .

Figure 1. Pareto optimal $\alpha_k$ (left and middle columns) and $\varphi$ (right column), with with $x=70$ , 0% interest, $\psi=0.1$ and $T=25$ . Each row provides the results for a different combination of $\gamma$ and $\delta$ . The left and middle columns display the optimal $\alpha_k$ for $N_0=100$ and $N_0=\infty$ , and within each figure, the red straight curves correspond to the $\alpha_k$ ’s for different values of $\nu\in(0,1)$ under the baseline systematic risk level with shades from dark to light corresponding to $\nu$ from 0 to 1. The dashed blue curves correspond to the high systematic risk scenario, with shades from dark to light corresponding to $\nu$ from 0 to 1. The right column displays the values of $\frac{\varphi}{\psi}$ with $\nu$ on the x-axis, where dots and squares are used for $N_0=100$ and $N_0=\infty$ , respectively, and gray and black are used for the baseline and high systematic risk levels, respectively.

The results indicate that the PO $\varphi$ typically takes values of either 0 (the policyholder’s optimal) or $\psi$ (the insurer’s optimal), with intermediary values emerging only when $\nu$ approaches 0.5, balancing the two objectives. This pattern is consistent with the individually optimal results, where the insurer and policyholder have opposing preferences for $\varphi$ .

The PO $\alpha_k$ typically lies between the individually optimal values for both parties. However, since in several cases and for some k’s the optimal $\alpha_k$ for the two parties are both equal to 1, the PO $\alpha_k$ is unambiguously equal to 1. The results in Figure 1, as well as in the other unreported combinations of input parameters, suggest that under the Lee–Carter model, the PO $\alpha_k$ has at most two inflection points over the course of the policy, with a shape determined predominantly by levels of systematic and diversifiable risk. For instance, for $N_0=100$ , a change from the baseline to the high systematic risk level appears to accentuate the inflection points. On the other hand, increasing portfolio size to $N_0=\infty$ appears to always lead to a U-shaped function. Note that under the CBD mortality model, the PO $\alpha_k$ can exhibit multiple inflection points in some cases.

To understand what drives the behavior of $\alpha_k$ over time, it is important to determine whether it is shaped by the preferences of policyholders or the insurer. Examining the evolution of $\alpha_k$ as a function of $\nu$ , it appears that as more weight is placed on the insurer’s objective, $\alpha_k$ tends to flatten toward 1. This indicates that the overall shape of the PO $\alpha_k$ is primarily driven by policyholders’ preferences rather than those of the insurer.

The shape of policyholder’s optimal $\alpha_k$ can be understood by analyzing how systematic and diversifiable risks interact across three different stages of the contract: later years (large k), mid-term (intermediate k), and early years (small k). First, toward the end of the contract, that is large values of k, survival probabilities are low, and the impact of risk is dampened, making it more acceptable for policyholders to take on greater risk-sharing. This pattern appears consistently across different levels of risk aversion and alternative mortality models (left and middle panels of Figure 1 for sufficiently high $\alpha_k$ ). Second, midway through the term, that is intermediate values of k, policyholders tend to avoid risk sharing at certain ages where survival probabilities’ variability is high and, unlike at older ages, the values of the probabilities are not high enough to offset the risk. This results in a local minimum in $\alpha_k$ , observed across most parameter choices. Last, at the beginning of the contract, that is small values of k, variability is driven more by portfolio size than by systematic risk. For small portfolios, higher volatility leads policyholders to avoid risk sharing ( $\alpha_k$ closer to 1 on the left panels of the figures), whereas, in large portfolios, lower diversifiable risk allows them to accept risk-sharing in exchange for reduced premiums ( $\alpha_k$ closer to 1 on the middle panels of the figures). Notably, increasing systematic risk has little effect on early-year risk-sharing decisions, while increasing portfolio size significantly shifts $\alpha_k$ from near zero to near one. Taken together, these effects explain why $\alpha_k$ follows a U-shape in large portfolios, with high risk sharing at the beginning and end and lower risk sharing in the middle. For smaller portfolios, the same mechanisms justify an increasing trend in $\alpha_k$ , potentially with an additional inflection point before the mid-term minimum in cases of high systematic risk.

The effect of the contract term T is analyzed separately as it highlights a pattern crucial to the remainder of the paper. To illustrate this pattern, it is helpful to first examine the individually optimal contracts, that is, for $\nu = 0$ and $\nu = 1$ . From the insurer’s viewpoint ( $\nu = 0$ ), the optimal $\alpha_k$ values are equal to 1 for all input parameter combinations tested, regardless of the contract term. From the policyholder’s viewpoint ( $\nu = 1$ ), Figure 2 shows the optimal $\alpha_k$ for a policyholder aged 70 years. The figure assumes 0% interest and $\psi = 0.1$ , but similar results hold for other input parameters, including for age 55.

Figure 2. Individually optimal $\alpha_k$ for the policyholder as a function k, with $x=70$ , 0% interest, and $\psi=0.1$ . Within each panel, the optimal $\alpha_k$ are reported for $T=25$ (pink), $T=20$ (red), $T=15$ (blue), $T=10$ (green) and $T=5$ (brown). The black curves are the combinations of all the $\alpha_k$ ’s over each interval. The solid curves correspond to the baseline systematic risk level, and the dashed curves correspond to the case with high systematic risk level. Each column gives the results for a different portfolio size, with $N_0=100$ on the left and $N_0=\infty$ on the right. Each row gives the result for a different combination of the risk aversion parameters $\gamma$ and $\delta$ .

To simplify the discussion, let $\alpha_k^{(T)}$ denote the policyholder’s share of risk for a contract term T. Each panel of Figure 2 displays $\alpha_k^{(25)}$ in pink, $\alpha_k^{(20)}$ in red, $\alpha_k^{(15)}$ in blue, $\alpha_k^{(10)}$ in green, and $\alpha_k^{(5)}$ in brown. The black curves represent a combination $\tilde{\alpha}_k$ , constructed such that:

\begin{align*}\tilde{\alpha}_k = \alpha_k^{(T_i)}, \quad \text{for } k = T_{i-1} + 1, \dots, T_i,\end{align*}

where $T_1, \dots, T_5$ are the terms 5, 10, …, 25, respectively. In Figure 2, the solid curves represent the baseline systematic risk level, and the dashed curves represent high systematic risk. Each column corresponds to a different portfolio size, and each row shows results for different combinations of the risk aversion parameters $\gamma$ and $\delta$ .

The main conclusion is that the contract term T has little effect on the optimal $\alpha_k$ . Specifically, $\alpha_k^{(10)}, ..., \alpha_k^{(25)}$ for $k = 1, \dots, 5$ can often all be approximated by $\alpha_k^{(5)}$ for $k = 1, \dots, 5$ . More generally:

\begin{align*}\alpha_k^{(T_i)} \approx \alpha_k^{(T_j)}, \quad \text{for } k = 1, \dots, T_j,\end{align*}

where $T_j \leq T_i$ . This approximation is less precise in the top-left panel (high systematic risk, high diversifiable risk, and low risk aversions). But even in these cases, it still provides a reasonable approximation of $\alpha$ ’s shape. The invariance of $\alpha_k$ to the time horizon is likely due to the definition of risk sharing, which is based on reserve deviations within each year. As a result, $\alpha_k$ at a given time depends primarily on the risk experience of that specific year, rather than the total contract length.

To further validate this pattern, but for more general PO designs, the regression of the form $\alpha^{(T_i)}= c\alpha^{(T_{j})}+\varepsilon$ is estimated, with $T_i\gt T_j$ , where $\alpha^{(T_i)}$ and $\alpha^{(T_{j})}$ are the PO policyholder’s share for all combinations of input parameters and $\nu$ . The regression is estimated for $T_i\in\{10,15,...,25\}$ where $T_{j}\in\{5,10,...,20\}$ .

Table 1 presents the estimated coefficient c for all combinations of $T_i$ and $T_j$ , with the adjusted R $^2$ in brackets. The results strongly support the adequacy of the approximation $\alpha_k^{(T_i)} \approx \alpha_k^{(T_j)}$ , as the estimated parameter c is close to 1 for all combinations of $T_i$ and $T_j$ , and the adjusted R $^2$ is close to 1. These results further highlight that the term T has a minimal effect on the shape of the PO $\alpha_k$ , allowing for effective approximations from shorter-term policies to longer-term ones. Note that similar results were found under the CBD model.

Table 1. Estimated parameter c (with adjusted $R^2$ in brackets) of the regression model $\alpha^{(T_i)}= c\alpha^{(T_{j})}+\varepsilon$ , with $T_i\gt T_j$ , where $\alpha^{(T_i)}$ and $\alpha^{(T_{j})}$ are the Pareto optimal shares for all combinations of input parameters (systematic risk level, risk aversions, initial loading, interest rate and age) and all values of $\nu\in(0,1)$ .

To sum up, this case study leads to three major conclusions regarding the PO risk sharing. First, the PO $\varphi$ is mostly equal to the policyholder’s optimal $\varphi$ for $\nu$ close to 1, or the insurer’s optimal $\varphi$ for $\nu$ close to 0, with intermediary values when $\nu$ is around 0.5. Second, the PO $\alpha_k$ can vary significantly between 0 and 1, and possibly with multiple inflection points, depending on the underlying parameters. Third, the results show that the term T of the policy does not affect the PO $\alpha_k$ ’s. Specifically, for $k\leq T_{j}\lt T_i$ , the PO $\alpha_k^{(T_i)}$ for a policy with term $T_i$ can reasonably be approximated by the PO $\alpha_k^{(T_j)}$ for a policy with term $T_j$ .

3.3. Closed-form heuristic for Pareto optimal risk sharing

This subsection develops a heuristic for determining the PO risk-sharing parameters $\underline{\alpha}$ and $\varphi$ in closed form, based on insights from the numerical analysis. The heuristic consists of a sequential approach, where shorter-term policies are used to approximate longer-term ones.

The motivation for this heuristic stems from a key numerical finding: the contract term T has little influence on the shape of the optimal $\underline{\alpha}$ . Specifically, the results indicate that for any given time k, the policyholder’s share of risk $\alpha_k$ in a long-term contract can be well approximated by the corresponding $\alpha_k$ in a shorter-term contract. Supported by regression analysis, this pattern holds across a wide range of input parameters, with the exception of certain extreme cases. The intuition behind this is that the risk-sharing mechanism is defined in terms of reserve deviations per year, meaning that yearly deviations contribute more than their compounded effect over the contract length.

Given these observations, a sequential optimization approach emerges as a natural way to approximate the PO solution. The heuristic begins by determining the PO values of $\varphi$ and $\alpha_1$ for a policy term of 1. Once these values are identified, they are then used to calculate $\alpha_2$ for a term of 2, with $\varphi$ and $\alpha_1$ held constant. This sequential process continues by finding $\alpha_k$ for a policy term of $k\leq T$ based on the values of $\varphi$ and $\alpha_1, \dots, \alpha_{k-1}$ from the preceding steps.

Using a superscript to highlight the dependence on the term, let $\mathcal{F}^{(1)}_{\nu}(\alpha_1,\varphi)$ be the objective function for a term 1, and $\mathcal{F}^{(k)}_{\nu}(\alpha_1,...,\alpha_{k-1},\varphi;\alpha_k)$ be the objective function for terms $2,..., T$ with $\alpha_1,...,\alpha_{k-1}$ and $\varphi$ held constant. The heuristic can be formulated as follows:

(3.4) \begin{eqnarray} \left(\alpha_1^*,\varphi^*\right) &=&\underset{\left(\alpha_1,\varphi\right)\in[0,1]\times[0,\psi]}{\arg\max} \mathcal{F}^{(1)}_{\nu}\left(\alpha_1,\varphi\right),\end{eqnarray}

(3.5) \begin{eqnarray} \alpha_k^*&=& \underset{\alpha_k\in[0,1]}{\arg\max} \quad \mathcal{F}^{(k)}_{\nu}\left(\alpha_1^*,...,\alpha_{k-1}^*,\varphi^*;\alpha_k\right), \qquad k=2,...,T.\end{eqnarray}

The following algorithmic pseudo-code summarizes the sequential optimization process to further clarify the stepwise approach:

The closed-form expression of the heuristic PO risk sharing is provided in the following proposition. The proof is in Appendix A.3.

Proposition 2. Consider the quantities $\Lambda_k^{P}=2\gamma b\sigma_{p,a}^{(k)}-\mu^{(k,k)}$ , $\Lambda_k^{I}=2\delta b\tilde{\sigma}_{p,a}^{(k)}+\mu^{(k,k)}$ , $\Lambda_k=\nu\Lambda_k^P+$ $(1-\nu)\Lambda_k^I$ , as well as

\begin{eqnarray*} \Gamma_{\nu}^{(\,j,k)} &=&\nu\gamma \sigma_p^{(\,j,k)} \alpha_j^*- (1-\nu) \delta \tilde{\sigma}_p^{(\,j,k)}(1-\alpha_j^*),\\ \Gamma_{\nu}^{(k)} &=&\nu\gamma \sigma_p^{(k,k)}- (1-\nu) \delta \tilde{\sigma}_p^{(k,k)}. \end{eqnarray*}

Consider also the functions $F(\alpha,\varphi)$ and $G(\alpha,\varphi)$ given by

\begin{eqnarray*} F(\alpha,\varphi)&=&b(1+\varphi)\left(\Lambda_1 - 2b(1+\varphi)\left(\nu \gamma \sigma_p^{(1,1)}\alpha - (1-\nu)\delta\tilde{\sigma}_p^{(1,1)}(1-\alpha)\right)\right),\\ G(\alpha,\varphi)&=&b\big((1-2\nu)a_{x:\overline{1}\kern-1.4pt\raise-1.2pt\hbox{$|$}} + \nu \alpha\Lambda_1^P - (1-\nu)(1-\alpha)\Lambda_1^I - 2b(1+\varphi)\mathcal{S}_\nu(\alpha)\big), \end{eqnarray*}

where $\mathcal{S}_\nu(\alpha)=\nu\gamma \alpha^2\sigma_p^{(1,1)} + (1-\nu)\delta (1-\alpha)^2\tilde{\sigma}_p^{(1,1)}$ .

Under the sequential heuristic PO risk sharing, $\alpha_1^*\in[0,1]$ and $\varphi^*\in[0,\psi]$ satisfy one of the following:

• • $\alpha_1^*=0$ and $\varphi^*=0$ , with $F(0,0)\lt0$ and $G(0,0)\lt0$ ;

• • $\alpha_1^*=0$ and $\varphi^*=\psi$ , with $F(0,\psi)\lt0$ and $G(0,\psi)\gt0$ ;

• • $\alpha_1^*=1$ and $\varphi^*=0$ , with $F(1,0)\gt0$ and $G(1,0)\lt0$ ;

• • $\alpha_1^*=1$ and $\varphi^*=\psi$ , with $F(1,\psi)\gt0$ and $G(1,\psi)\gt0$ ;

• • $\alpha^*_1=\frac{\frac{\Lambda_1}{2b} + (1-\nu)\delta \tilde{\sigma}_p^{(1,1)}}{\nu\gamma\sigma_p^{(1,1)} + (1-\nu)\delta \tilde{\sigma}_p^{(1,1)}}$ and $\varphi^*=0$ , with $G(\alpha_1^*,0)\lt0$ ;

• • $\alpha^*_1=\frac{\frac{\Lambda_1}{2b(1+\psi)} + (1-\nu)\delta \tilde{\sigma}_p^{(1,1)}}{\nu\gamma\sigma_p^{(1,1)} + (1-\nu)\delta \tilde{\sigma}_p^{(1,1)}}$ and $\varphi^*=\psi$ , with $G(\alpha_1^*,\psi)\gt0$ ;

• • $\alpha_1^*=0$ and $\varphi^*=\frac{a_{x:\overline{1}\kern-1.4pt\raise-2.5pt\hbox{$|$}}(1-2\nu) - (1-\nu)\Lambda_1^I}{2b(1-\nu)\delta\tilde{\sigma}_p^{(1,1)}}-1$ , with $F(0,\varphi^*)\lt0$ ;

• • $\alpha_1^*=1$ and $\varphi^*=\frac{a_{x:\overline{1}\kern-1.4pt\raise-2.5pt\hbox{$|$}}(1-2\nu) + \nu\Lambda_1^P}{2b\nu\gamma\sigma_p^{(1,1)}}-1$ , with $F(1,\varphi^*)\gt0$ ;

• • $\alpha_1^*=\frac{\Lambda_1(1-\nu)\delta \tilde{\sigma}_p^{(1,1)}+(1-\nu)\delta\tilde{\sigma}_p^{(1,1)}\big(a_{x:\overline{1}\kern-1.4pt\raise-2.5pt\hbox{$|$}}(1-2\nu)-(1-\nu)\Lambda_1^I\big)}{\left(\nu\gamma\sigma_p^{(1,1)}+(1-\nu)\delta\tilde{\sigma}_p^{(1,1)}\right)\big(a_{x:\overline{1}\kern-1.4pt\raise-2.5pt\hbox{$|$}}(1-2\nu)-(1-\nu)\Lambda_1^I\big)+\Lambda_1(1-\nu)\delta \tilde{\sigma}_p^{(1,1)}}$ and $\varphi^* = \frac{a_{x:\overline{1}\kern-1.4pt\raise-2.5pt\hbox{$|$}}(1-2\nu)-(1-\nu)\Lambda_1^I}{2b(1-\nu)\delta \tilde{\sigma}_p^{(1,1)}(1-\alpha_1^*)}-1$ .

Furthermore, for $k=2,...,T$ , and given $\varphi^*$ and $\alpha_1^*,...,\alpha_{k-1}^*$ , the optimal $\alpha_k^*\in[0,1]$ satisfies one of the following:

• • $\alpha_k^*=0$ , with $\Lambda_k+2b(1-\nu)(1+\varphi^*) \delta \tilde{\sigma}_p^{(k,k)} - 2b(1+\varphi^*)\sum_{j=1}^{k-1}\Gamma_{\nu}^{(\,j,k)}\lt0$ ;

• • $\alpha_k^*=1$ , with $\Lambda_k+2b(1-\nu)(1+\varphi^*) \delta \tilde{\sigma}_p^{(k,k)} - 2b(1+\varphi^*)\sum_{j=1}^{k-1}\Gamma_{\nu}^{(\,j,k)}- 2b(1+\varphi^*)\Gamma_{\nu}^{(k)}\gt0$ ;

• • $ \alpha_k^*=\frac{1}{\Gamma_{\nu}^{(k)}}\left(\frac{\Lambda_k}{2b(1+\varphi^*)}+(1-\nu)\delta \tilde{\sigma}_p^{(k,k)} - \sum_{j=1}^{k-1}\Gamma_{\nu}^{(\,j,k)}\right).$

Proposition 2 simplifies the optimization problem by solving for $\alpha_1$ and $\varphi$ for a single-term policy and then iteratively determining subsequent $\alpha_k$ values for longer terms. This sequential approach substantially reduces the computational complexity of the full optimization problem. Moreover, it provides a more manageable framework to account for the constraints, where constraints are handled at each stage of the policy term.

Figure 3 displays the relative difference in percentage between the optimal objective function values obtained from the numerical optimization and those from the heuristic. Specifically, for $\left(\underline{\hat{\alpha}},\hat{\varphi}\right)$ the parameter values estimated using the numerical approach, and $\left(\underline{\alpha}^*,\varphi^*\right)$ those obtained from the heuristic in Proposition 2, let $\hat{\mathcal{F}}_\nu=\mathcal{F}_\nu\left(\underline{\hat{\alpha}},\hat{\varphi}\right)$ and $\mathcal{F}_\nu^*=\mathcal{F}_\nu\left(\underline{\alpha}^*,\varphi^*\right)$ . The figure displays $100\times\frac{\hat{\mathcal{F}}_\nu-\mathcal{F}_\nu^*}{|\hat{\mathcal{F}}_\nu|}$ as a function of $\nu$ , with 0% interest, $\psi=0.1$ and a term $T=25$ under different combinations of age, risk aversion parameters, portfolio size, and systematic risk levels.

Figure 3. Relative difference $100\times\frac{\hat{\mathcal{F}}_\nu-\mathcal{F}_\nu^*}{|\hat{\mathcal{F}}_\nu|}$ between the optimal objective functions obtained from the numerical optimization $\hat{\mathcal{F}}_\nu$ and that obtained from the heuristic $\mathcal{F}_\nu^*$ , with 0% interest, $\psi=0.1$ and a term $T=25$ . Within each panel, the relative difference is given in function of $\nu$ for a policyholder aged $x=55$ in black and $x=70$ in red. The solid curves correspond to the baseline systematic risk level, and the dashed curves correspond to the case with high systematic risk level. Each column gives the results for a different portfolio size, with $N_0=100$ on the left and $N_0=\infty$ on the right. Each row gives the result for a different combination of the risk aversion parameters $\gamma$ and $\delta$ .

Consistent with earlier observations that the policy term has minimal impact, the figure demonstrates that the heuristic provides a very close approximation to the numerical solution, with relative differences below 0.1% for nearly all combinations of input parameters. Unreported results for a 5% interest rate show similarly small differences. Increasing $\varphi$ ’s upper bound $\psi$ from 0.1 to 0.3 results in slightly larger, but still very low, relative differences, with a maximum of about 1% under high systematic risk and low risk aversion parameters $\gamma$ and $\delta$ . In some cases, however, where the objective function was very close to 0, the relative difference was disproportionately high, with one instance slightly exceeding 100%. Despite this, the heuristic remains effective, as these discrepancies occur only when the objective function is near zero and thus less impactful.

Finally, Figure 4 displays $\frac{1}{T}\sum_{k=1}^{T}|\hat{\alpha}_k-\alpha_k^*|$ , that is, the absolute difference averaged over the term between $\hat{\alpha}_k$ obtained numerically and $\alpha^*_k$ obtained from the heuristic. The figure shows small differences between $\hat{\alpha}_k$ and $\alpha_k^*$ , which tend to increase when $\nu$ is above 0.5. Overall, along with the results in Figure 3, the differences between $\hat{\alpha}_k$ and $\alpha_k^*$ displayed in Figure 4 show that the heuristic provides a good approximation for the PO risk sharing.

Figure 4. Mean absolute difference $\frac{1}{T}\sum_{k=1}^{T}|\hat{\alpha}_k-\alpha_k^*|$ between the optimal $\alpha_k$ ’s averaged over the term, where $\hat{\alpha}_k$ are those obtained from the numerical optimization and $\alpha^*_k$ are those obtained from the heuristic, with 0% interest, $\psi=0.1$ and a term $T=25$ . Within each panel, the absolute difference is given in function of $\nu$ for a policyholder aged $x=55$ in black and $x=70$ in red. The solid curves correspond to the baseline systematic risk level, and the dashed curves correspond to the case with high systematic risk level. Each column gives the results for a different portfolio size, with $N_0=100$ on the left and $N_0=\infty$ on the right. Each row gives the result for a different combination of the risk aversion parameters $\gamma$ and  $\delta$ .

3.4. Comparison of constant and heuristic time-varying risk sharing

Before concluding this section, it is natural to ask whether the added complexity of time-varying risk sharing provides meaningful benefits compared to the simpler constant risk sharing, which is a special case of the time-varying approach. This question is particularly relevant for the heuristic risk-sharing rule, as it approximates the PO time-varying strategy but does not necessarily guarantee better outcomes than the constant case. The key issue is which party benefits from time variation in risk sharing under the heuristic rule and the extent of these gains or losses.

To address this, individual preferences under a constant risk-sharing strategy are compared with those under the heuristic time-varying risk-sharing approach. Specifically, the percentage gains in preference from adopting the heuristic time-varying strategy are given by

\begin{align*}100\times\frac{\rho^{pol}(\underline{\alpha}^{*},\varphi^*) - \rho^{pol}(\alpha^c,\varphi^c)}{|\rho^{pol}(\alpha^c,\varphi^c)|} \quad \text{ and } \quad 100 \times \frac{\rho^{ins}(\underline{\alpha}^{*},\varphi^*) - \rho^{ins}(\alpha^c,\varphi^c)}{|\rho^{ins}(\alpha^c,\varphi^c)|},\end{align*}

for the policyholder and the insurer, respectively, where $\underline{\alpha}^{*}$ and $\varphi^*$ represent the heuristic time-varying PO risk-sharing coefficients and initial loading, while $\alpha^c$ and $\varphi^c$ correspond to their counterparts under constant PO risk sharing.

Figure 5 provides the results for different parameter settings for the policyholder (left panels) and the insurer (right panels). Each row represents a different level of risk aversion, and each panel shows how the gains evolve as a function of $\nu$ , under different combinations of systematic risk and portfolio size.

Figure 5. Preference gains, calculated as $100\times\frac{\rho^{pol}(\underline{\alpha}^{*},\varphi^*) - \rho^{pol}(\alpha^c,\varphi^c)}{|\rho^{pol}(\alpha^c,\varphi^c)|}$ and $100 \times \frac{\rho^{ins}(\underline{\alpha}^{*},\varphi^*) - \rho^{ins}(\alpha^c,\varphi^c)}{|\rho^{ins}(\alpha^c,\varphi^c)|}$ for the policyholder (left panels) and insurer (right panels), respectively, when moving from constant to heuristic time-varying risk-sharing. Results are obtained for age $x=70$ , 0% interest, $\psi=0.1$ and a term $T=25$ . Within each panel, the relative difference is given in function of $\nu$ . The solid curves correspond to the baseline systematic risk level, and the dashed curves correspond to the case with high systematic risk level. Black curves correspond to the case of high diversifiable risk ( $N_0=100$ ), and red ones correspond to the case of low diversifiable risk ( $N_0=\infty$ ). Each row gives the result for a different combination of the risk aversion parameters $\gamma$ and $\delta$ .

The results indicate that for policyholders, the shift to time-varying risk sharing does not lead to a substantial improvement or deterioration in preferences across all tested cases, suggesting that the added complexity does not significantly affect their expected utility. In contrast, for insurers, the effect is consistently positive and often substantial, with preference gains of about 20% in many cases, and exceeding 100% in some instances. This highlights that the heuristic time-varying risk sharing provides a significant advantage to the insurer compared to the constant risk sharing, without imposing any meaningful cost on the policyholder.