2. The modeling framework

Let $(\Omega, \mathcal{G}, (\mathcal{G}_t )_t, \mathbb{P})$ denote a filtered probability space where $\mathbb{P}$ is the historical probability measure. We define $\tau_x$ and $\tau_y$ as future lifetimes for annuitants aged x and y, respectively. These lifetimes have corresponding stochastic intensities $(\mu_x(t))_{t\geq 0}$ and $(\mu_y(t))_{t\geq 0}$ . Referring to the annuitant aged x as x and aged y as y, and inspired by Dahl (Reference Dahl2004) and Biffis (Reference Biffis2005), $\tau_x$ marks the initial jump-time of a $\mathcal{G}$ -counting process $(N_t^x)_{t\geq 0}$ , indicating if individual x has passed away $(N_t^x \neq 0)$ by time t. Similarly, $\tau_y$ indicates the same for individual y.

Drawing from Deelstra et al. (Reference Deelstra, Grasselli and Van Weverberg2016), $\mathcal{R}_t$ and $\mathcal{M}_t$ are filtrations produced by the interest rate and joint mortality intensity processes, respectively. Let $\mathcal{F}_t \, :\!= \, \mathcal{R}_t \vee \mathcal{M}_t$ be the minimal $\sigma$ -algebra that contains $\mathcal{R}_t \cup \mathcal{M}_t$ . The filtration $(\mathcal{G}_t)_{t\geq 0}$ embodies the progressive information available, capturing the evolution of both state variables and whether the annuitants x or y have died at time t. Hence, we set $\mathcal{G}_t \, :\!= \, \mathcal{F}_t \vee \mathcal{H}_t$ , where $\mathcal{H}_t\, :\!= \, \sigma(\{\{\tau_x \leq s\},\, \{\tau_y \leq s\} ,\,0 \leq s \leq t\})$ is the minimum filtration wherein $\tau_x$ and $\tau_y$ are established stopping times.

In line with Deelstra et al. (Reference Deelstra, Grasselli and Van Weverberg2016), we posit that $(N_t^x,N_t^y)_{t\geq 0}$ operates as a two-dimensional Cox process, driven by an $\mathcal{F}$ sub-filtration of $\mathcal{G}$ . This implies the following survival probability on the set $\lbrace \tau_x \gt t, \tau_y \gt t \rbrace $ :

(1) \begin{align}\mathbb{Q}(\tau_x\gt T, \tau_y\gt T|\mathcal{G}_t)=\mathbb{E}^{\mathbb{Q}}\left[\left. e^{-\int_t^T\mu_x(s) + \mu_y(s)ds} \right\vert \mathcal{G}_t\right]\!, \end{align}

where $\mathbb{Q}$ is a risk-neutral pricing measure equivalent to $\mathbb{P}$ and $\mathbb{E}^{\mathbb{Q}}\left[\cdot|\mathcal{G}_t\right]$ denotes the conditional expectation under $\mathbb{Q}$ given $\mathcal{G}_t$ . Regarding the choice of $\mathbb{Q}$ , we follow Deelstra et al. (Reference Deelstra, Grasselli and Van Weverberg2016) and references therein and assume that it can be selected on the basis of available market data. For literature addressing the non-diversifiable aspect of mortality risk, which results in an incomplete market, refer to studies such as Dahl and Mø ller (Reference Dahl and Møller2006), Bayraktar et al. (Reference Bayraktar, Milevsky, Promislow and Young2009), Huang et al. (Reference Huang, Tsai, Yang and Cheng2014), and Ceci et al. (Reference Ceci, Colaneri and Cretarola2015).

For completeness, we restate some well-known equations in the literature. We note that since $(N_t^x, N_t^y)$ is a two-dimensional Cox process, then it follows that, conditional on a particular trajectory $(\mu_x(t),\mu_y(t))_{t \geq 0}$ , the process $(N_t^x,N_t^y)_{t \geq 0}$ is a two-dimensional independent Poisson-inhomogeneous process with parameter $(\!\int_0^t\mu_x(s)ds)_{t \geq 0}$ for $(N_t^x)_{t \geq 0}$ and $(\!\int_0^t\mu_y(s)ds)_{t \geq 0}$ for $(N_t^y)_{t \geq 0}$ . As such, for all $t_1 \geq t \geq 0 $ and $t_2 \geq t \geq 0 $ and non-negative integers $k_1$ , $k_2$ , the following equality holds

\begin{align*}I(k_1,k_2) \, :\!= \,& \,\mathbb{Q}( (N_{t_1}^x - N_t^x = k_1) (N_{t_2}^y- N_t^y = k_2)| \mathcal{H}_{t} \vee \mathcal{F}_{\infty})\\[3pt] =& \,\mathbb{Q}( (N_{t_1}^x - N_t^x = k_1) | \mathcal{H}_{t} \vee \mathcal{F}_{\infty}) \mathbb{Q}( (N_{t_2}^y- N_t^y = k_2)| \mathcal{H}_{t} \vee \mathcal{F}_{\infty})\\[3pt] =& \,\frac{\left(\int_t^{t_1}\mu_x(s)ds\right)^{k_1} }{k_1!}e^{-\int_t^{t_1}\mu_x(s)ds}\frac{\left(\int_t^{t_2}\mu_y(s)ds\right)^{k_2} }{k_2!}e^{-\int_t^{t_2}\mu_y(s)ds},\end{align*}

where we used the conditional independence between $(N_t^{x})_{t\geq 0}$ and $(N_t^{y})_{t\geq 0}$ . Taking $k_1=0$ and $k_2=0$ and integrating the intensities lead to (1).

Still following the literature, see for example Dahl (Reference Dahl2004, Section 3.4) and Biffis (Reference Biffis2005, Equation (10)), the $\mathcal{G}_t$ -conditional density of $(\tau_x,\tau_y)$ on the set $\lbrace \tau_x\gt t, \tau_y\gt t\rbrace$ , denoted $f_t(t_1,t_2)$ , is given by $f_t(t_1,t_2) = \partial^2_{t_1t_2}\mathbb{Q}( (\tau_x\gt t_1)(\tau_y\gt t_2) | \mathcal{G}_t)$ . Note that nothing so far was said on the intensities $(\mu_x(t),\mu_y(t))_{t \geq 0}$ apart from the fact that they have to be positive. If they are dependent, then there is a dependency between the mortality intensities even if conditionally on the intensities $(N_t^x)_{t \geq 0}$ and $(N_t^y)_{t \geq 0}$ are independent. It is through the intensities that we can incorporate a dependency between the mortality events $\tau_x$ and $\tau_y$ .

Assuming $0 \leq t \leq T$ , then the pricing at time t of an insurance contingent claim payable at T paying 1 on the survival of annuitants aged x and y at time 0 is given on the set $\lbrace \tau_x\gt t, \tau_y\gt t \rbrace $ by

(2) \begin{align}\mathrm{SB}(t, T)\, :\!= \,\mathbb{E}^{\mathbb{Q}}\left[\left. e^{-\int_t^T r(s)ds}\mathbf{1}_{\lbrace \tau_x\gt T \rbrace}\mathbf{1}_{\lbrace \tau_y\gt T \rbrace}\right\vert \mathcal{G}_t \right], \end{align}

where $(r(t))_{t\geq 0}$ is the risk-free rate which is adapted to the filtration $\left(\mathcal{F}_t\right)_{t\geq 0}$ . This instrument is commonly called the joint-survival zero-coupon bond. At this stage, it is required to specify whether the risk-free rate $(r(t))_{t\geq 0}$ depends on the mortality intensities $(\mu_x(t),\mu_y(t))_{t \geq 0}$ . We follow the literature and assume that interest rates and the mortality risks are independent from each other, see for example Biffis (Reference Biffis2005).Footnote ¹

Furthermore, using (1) as well as the fact that the conditioning on $\mathcal{G}_t$ can be reduced to that on $\mathcal{F}_t$ as shown in Biffis (Reference Biffis2005, Appendix C), we obtain

(3) \begin{align}\mathrm{SB}(t, T)=\mathbf{1}_{\lbrace \tau_x\gt t,\tau_y\gt t \rbrace}P(t,T)\mathbb{E}^{\mathbb{Q}}\left[\left. e^{-\int_t^T \mu_x(s) + \mu_y(s)ds}\right\vert \mathcal{F}_t \right], \end{align}

with P(t,T) the time t-value of a zero-coupon bond with maturity T. Hereafter, $\mathbb{E}^{\mathbb{Q}}\left[\cdot| \mathcal{F}_t \right]$ will be denoted $\mathbb{E}_t^{\mathbb{Q}}\left[\cdot\right]$ so that the joint-survival (zero-coupon) bond, denoted by $\mathrm{SB}(t,T)$ , is eventually given on the set $\lbrace \tau_x\gt t, \tau_y\gt t \rbrace $ by

(4) \begin{align}\mathrm{SB}(t,T) &= P(t,T)\mathbb{E}_t^{\mathbb{Q}}\left[e^{-\int_t^T \mu_x(s)+ \mu_y(s) ds}\right]= P(t,T)\mathrm{SB}_0(t,T), \end{align}

where we use $\mathrm{SB}_0(\cdot, \cdot)$ to denote the joint-survival bond that contains no interest rate risk. The equivalent actuarial notation is given by ${}_Tp_{xy}$ ; however, due to the similarities of our work with Biffis (Reference Biffis2005) and Deelstra et al. (Reference Deelstra, Grasselli and Van Weverberg2016), we instead follow their notation.

Using the framework introduced so far, the joint-life annuity is therefore the sum of joint-survival bonds, that is, a product that pays at certain prespecified dates an amount (possibly date dependent), conditional on the survival of both annuitants at those dates. Let $T_1,\ldots,T_N$ be a set of yearly spaced dates such that $T_N=\min\!(x^*-x-1,y^*-y-1)$ with $x^*$ and $y^*$ being the maximum attainable age for the annuitants x and y, respectively. Consider the product that pays at times $T_i$ the amount 1 conditional on the survival of both annuitants, the value of the product at time t is given on the set $\lbrace \tau_x\gt t, \tau_y\gt t \rbrace $ by

(5) \begin{align}\sum_{i=1}^N \mathrm{SB}(t,T_i), \end{align}

with $\mathrm{SB}(t,T)$ given by (4).

Another product of interest is the guaranteed joint-survival annuity option, that is, an option on (5). We follow Deelstra et al. (Reference Deelstra, Grasselli and Van Weverberg2016, Section 3.2), who analyze the guaranteed survival annuity option (for a single annuitant) but consider the case of multiple annuitants. Suppose a guaranteed joint-survival annuity option with maturity T and consider a set of yearly spaced dates $T_1,\ldots, T_N$ such that $T_1=T + 1$ (where the units are years) and $T_N=\min\!(x^*-(x+T)-1,y^*-(y+T)-1)$ with $x^*$ and $y^*$ again being the maximum attainable age for annuitants x and y, respectively. Then we denote the value of a guaranteed joint-survival annuity option at time t, exercisable at T, on a joint-survival annuity paying at times $T_1,...,T_N$ as $C(t,T,T_N)$ . Thus at $t=T$ , the payoff of the option is given on the set $\lbrace \tau_x\gt T, \tau_y\gt T \rbrace $  by

(6) \begin{align}C(t=T,T,T_N) &\, :\!= \, \max\left(g \sum_{i=1}^N \mathrm{SB}(T,T_i) ,1 \right) =1 + g\,\bar{\!C}(T,T,T_N), \end{align}

where g is the fixed rate called the guaranteed joint-survival annuity rate and $(z)_+=\max\!(z,0)$ and $\,\bar{\!C}(T,T,T_N) = \left(\sum_{i=1}^N \mathrm{SB}(T,T_i) -1/g\right)_+$ . Therefore, the guaranteed joint-survival annuity option value at time t is given on the set $\lbrace \tau_x\gt t, \tau_y\gt t \rbrace $ by

\begin{align*}C(t,T,T_N) &=P(t,T)\mathrm{SB}_0(t,T) + g \mathbb{E}_t^{\mathbb{Q}}\left[ e^{-\int_{t}^{T}(r(s) + \mu_x(s)+\mu_y(s))ds} \,\bar{\!C}(T,T,T_N) \right],\end{align*}

where $\,\bar{\!C}(t,T,T_N) $ is the expectation on the right-hand side above. It is worth having a closer look at the expression of $\,\bar{\!C}(t,T,T_N)$ as it highlights some modeling issues. Using (4), we obtain on the set $\lbrace \tau_x\gt t, \tau_y\gt t \rbrace $ the equality

(7) \begin{align}\,\bar{\!C}(t,T,T_N) = \mathbb{E}_t^{\mathbb{Q}}\left[ e^{-\int_{t}^{T}(r(s) + \mu_x(s)+\mu_y(s))ds}\left( \sum_{i=1}^NP(T,T_i)\mathrm{SB}_0(T,T_i)-1/g \right)_+ \right], \end{align}

where we can clearly see that even if the interest rate is independent with respect to the mortality intensities, its distribution impacts the option price. Given the difficulty of treating stochastic interest rates, independent or dependent of the mortality intensities, we take the interest rate to be deterministic.

Let us stress the fact that all of the equations stated above are model independent and are not contingent on using any specific model or process.

As seen from (7), the pricing of a guaranteed (joint) annuity option is a difficult task, as it requires the law of the sum of joint-survival (zero-coupon) bonds. If the joint-survival bond is an exponential function of the state variable, then this amounts to requiring the law of the sum of exponentials of random variables that is in general unknown. This difficulty is well explained in Biffis and Millossovich (Reference Biffis and Millossovich2006) and is one of the main difficulties we need to overcome.

3. The linear-rational Wishart model

To build the joint mortality model, we use the Wishart process defined in Bru (Reference Bru1991) and introduce its main properties. In the following section, we provide a concise review of the essential techniques and findings required to construct our model. This presentation is prompted by the relatively novel application of the Wishart process in mortality modeling. Given a filtered probability space $(\Omega,\mathcal{F}, (\mathcal{F}_t)_t, \mathbb{P})$ , we denote by $\mathbb{E}\left[\cdot \right]$ (resp. $\mathbb{E}_t\left[ \cdot \right]\, :\!= \,\mathbb{E}\left[\,\cdot\,|\mathcal{F}_t \right]$ ) the expectation (resp. conditional expectation) under the historical probability measure $\mathbb{P}$ . The Wishart process satisfies the matrix SDE

(8) \begin{equation}dv_t = (\omega + m v_t + v_t m^\top )dt + \sqrt{v_t}dw_t \sigma + \sigma^\top dw_t^\top \sqrt{v_t}\,, \end{equation}

where $v_t$ is an $n\times n$ matrix that belongs to the set of positive-definite matrices denoted $\mathbb{S}_n^{++}$ , $m,\sigma$ belong to the set of $n\times n$ real matrices denoted $\mathsf{M}(n)$ , $\lbrace w_t;\, t\geq 0 \rbrace$ is a matrix Brownian motion of dimension $n\times n$ (i.e., a matrix of $n^2$ independent scalar Brownian motions) under the probability measure $\mathbb{P}$ , and $\cdot^\top$ stands for the matrix transposition. The matrix $\omega \in \mathbb{S}_n^{++}$ satisfies certain constraints involving $\sigma^{\top}\sigma$ to ensure the positive-definiteness of the matrix process $v_t$ . The quantity $\sqrt{v_t}$ is well defined since $v_t \in \mathbb{S}_n^{++}$ . The matrix m is such that $\lbrace \Re(\lambda_i^{m})\lt 0;\, i=1,\ldots,n \rbrace $ where $\lambda_i^{m}\in \mathsf{Spec}(m)$ for $i=1,\ldots,n$ and $\mathsf{Spec}(m)$ is the spectrum of the matrix m, while $\Re(\,\cdot\,)$ denotes the real component. The matrix $\sigma$ belongs to $\mathsf{GL}_n(\mathbb{R})$ the general linear group over $\mathbb{R}$ (i.e., the set of real invertible matrices). Due to the invariance of the law of the Brownian motion with respect to rotations and the polar decomposition of $\sigma$ , we can assume that $\sigma \in \mathbb{S}_n^{++}$ . We denote $e_{ij}$ to be the basis of $\mathsf{M}(n)$ , that is, it is the $n\times n$ matrix with 1 in the (i, j) place and 0 elsewhere. Lastly $I_n$ stands for the $n\times n$ identity matrix, $0_n$ is the $n\times n$ null matrix, while $\text{diag}(z)$ with $z\in \mathbb{R}^n$ is the $n\times n$ matrix with z on its diagonal. The infinitesimal generator of the Wishart process is given by Bru (Reference Bru1991):

(9) \begin{align}\mathcal{L} = \text{tr}[(\omega + m v + vm^\top)\, D + 2 v D \sigma^2 D]\,, \end{align}

where $\text{tr}[\,\!\cdot\!\,]$ is the trace of a matrix and D is the $n\times n$ matrix operator $D_{ij}\, :\!= \,\partial_{v_{ij}}$ .

The Wishart process was initially defined and analyzed in Bru (Reference Bru1991) under the assumption that

(10) \begin{equation} \omega=\beta \sigma^2 \end{equation}

with $\beta \in \mathbb{R}_+$ such that $\beta \geq n+1$ to ensure that $v_t \in \mathbb{S}_n^{++}$ . Hereafter, this specification will be referred to as the Bru case. It was later extended in Mayerhofer et al. (Reference Mayerhofer, Pfaffel and Stelzer2011) (see also Cuchiero et al., Reference Cuchiero, Filipović, Mayerhofer and Teichmann2011) to the case $\omega \in \mathbb{S}_n^{++}$ and proved that if $ \omega - \beta \sigma^2 \in \mathbb{S}_n^{++}$ , with $\beta\geq n+1$ then $v_t\in \mathbb{S}_n^{++}$ .

In our work, we deliberately restrict ourselves to the Bru case (i.e., $\omega=\beta \sigma^2$ ). The advantage of the Bru case is that moment-generating function can be explicitly integrated making it fully explicit. This allows us to derive analytic formulae surrounding the guaranteed joint-survival annuity option and the individual mortalities in contrast to the work of Da Fonseca (Reference Da Fonseca2024) which is not in the Bru case. The below result is well known. We refer to Bru (Reference Bru1991) or Grasselli and Tebaldi (Reference Grasselli and Tebaldi2008) (see also Alfonsi, Reference Alfonsi2015, Chapter 5).Footnote ²

Proposition 3.1. Suppose that $\omega$ in (8) is such that $\omega=\beta \sigma^2$ with $\beta \in \mathbb{R}$ and $\beta \geq n+1$ , then the moment-generating function of $v_t$ for $\theta_1\in \mathbb{S}_n$ is given by

(11) \begin{align} \Phi(t,\theta_1,v_0) &=\mathbb{E}\left[ \exp\left( \text{tr}[\theta_1 v_t] \right)\right]=\frac{\text{etr}\left(\frac{\vartheta_t^\top}{2}(2\varsigma_t \theta_1)(I_n-2\varsigma_t\theta_1)^{-1} \right)}{\det(I_n -2\varsigma_t\theta_1)^{\beta/2}}, \end{align}

with $\varsigma_t \, :\!= \,\int_0^te^{(t-s)m}\sigma^2 e^{(t-s)m^\top}ds$ , which is given by $\text{vec}(\varsigma_t)=\mathsf{A}^{-1}(e^{\mathsf{A}t}-I_{n^2})\text{vec}(\sigma^2)$ , $\vartheta_t \, :\!= \, \varsigma_t^{-1}e^{mt}v_0e^{m^\top t}$ , $\text{etr}(\,\! \cdot\!\,)=\exp(\text{tr}[\,\!\cdot\!\,])$ and $\det(\! \cdot\!)$ the determinant of a matrix. Here, $\text{vec}(\!\cdot\!)$ denotes the vec operator that transforms a $n\times n$ matrix into a $n^2$ vector by stacking the columns, $\mathsf{A}\, :\!= \,I_{n}\otimes m+m\otimes I_{n}$ , and $\otimes$ denotes the Kronecker product.

The Wishart process is affine, that is, the moment-generating function is exponentially affine in the state variable. Equation (11) is the moment-generating function of a non-central Wishart distribution, hence the name of the process, and the corresponding density is given in Gupta and Nagar Reference Gupta and Nagar2000, Equation 3.5.1 p. 114). Note that even if the density is known, computing an expectation requires to perform a $n(n+1)/2$ -dimensional numerical integration.

The results developed in this work rely on the following result, which arises from solving Lyapunov equations after rewriting the dynamics of the Wishart process. This fact is well known in the matrix Riccati literature, see for example Abou-Kandil et al. (Reference Abou-Kandil, Freiling, Ionescu and Jank2003).

Lemma 3.3. Let $v_t$ be a Wishart process described by (8). If $u_0 \in \mathsf{M}(n)$ , then there exist a function $b_0(t)$ and a matrix function $a_0(t) \in \mathsf{M}(n)$ such that

(12) \begin{align}\mathbb{E}\left[\text{tr}[u_0v_{t}]\right]&=\text{tr}[a_0(t)v_{0}] + b_0(t), \end{align}

with

(13) \begin{align} \text{vec}\left(a_0(t)^{\top}\right) = e^{\mathsf{A}^{\top}t}\text{vec}(u_0^\top), \quad b_0(t) = \text{vec}(u_0^\top)^\top \mathsf{A}^{-1}\left(e^{\mathsf{A}t}-I_{n^2} \right)\mathsf{b}.\end{align}

Here, $\mathsf{b}\, :\!= \,\text{vec}\left(\omega\right)$ , $\mathsf{A}\, :\!= \,I_{n}\otimes m+m\otimes I_{n}$ , and $\otimes$ denotes the Kronecker product.

The following proposition is well known in the literature and is useful for obtaining the quadratic variation of linear functions of the Wishart process, which in turn will be utilized to derive the covariation (or instantaneous correlation) between mortality intensities. A proof of the proposition appears in p. 214 of Deelstra et al. (Reference Deelstra, Grasselli and Van Weverberg2016) and reads as follows.

Proposition 3.4. Given $h_1, h_2\in \mathsf{M}(n)$ the quadratic variation between $(\text{tr}[h_1v_t])_{t \geq 0 }$ and $(\text{tr}[h_2v_t])_{t \geq 0 }$ is

\begin{align*}d\langle \text{tr}[h_1v_.],\text{tr}[h_2v_.]\rangle_t&=\text{tr}[(h_1+h_1^\top) v_t(h_2+h_2^\top)\sigma^2] dt\,.\end{align*}

When $n=2$ , using Proposition 3.4, one can in particular show that the quadratic covariation between $v_{11,t}$ and $v_{22,t}$ is given by

(14) \begin{align}d\langle v_{11,.},v_{22,.}\rangle_t &= 4v_{12,t}(\sigma^2)_{12}dt ,\end{align}

where $(\sigma^2)_{ij}$ is the (i, j) element of the matrix $\sigma^2$ . The identity in (14) is specific to the two-dimensional case $n=2$ . Note that $(\sigma^2)_{12}=0$ if $\sigma_{12}=0$ , which means the dependence between $v_{11,t}$ and $v_{22,t}$ is primarily dependent on $\sigma_{12}$ (for $n=2$ ).

To obtain the derivative of the moment-generating function with respect to a parameter, the next proposition provides the details in the Bru case. This result will be particularly helpful for obtaining the densities of mortality intensities, as well as risk measures such as expected shortfall and value at risk for joint-survival annuities.

Proposition 3.5. Suppose that $\omega$ in (8) is such that $\omega=\beta \sigma^2$ with $\beta \in \mathbb{R}$ and $\beta \geq n+1$ and consider the moment-generating function of $v_t$ given by (11). For $\theta_1\in \mathbb{S}_n$ and $\theta_2\in \mathbb{S}_n$ consider the function for $\nu \in \mathbb{R}$ , $ \nu \to \Phi(t,\theta_1 + \nu \theta_2,v_0) = \mathbb{E}\left[\text{etr}\left( (\theta_1+\nu \theta_2) v_t \right)\right]$ then its derivative with respect to $\nu$ is given by

(15) \begin{align}\Phi_\nu(t,\theta_1,\theta_2,\nu, v_0) &\, :\!= \,\frac{d\Phi(t,\theta_1+\nu \theta_2, v_0)}{d\nu} =(g_1+g_2+g_3)\Phi(t,\theta_1+\nu \theta_2, v_0), \end{align}

with

\begin{align*}g_1(t,\theta_1,\theta_2,\nu) &= \text{tr}\left[\frac{\vartheta_t^\top}{2} 2\varsigma_t\theta_2(I_n -2\varsigma_t(\theta_1+\nu \theta_2))^{-1} \right], \\[3pt]g_2(t,\theta_1,\theta_2,\nu) &= \text{tr}\left[\frac{\vartheta_t^\top}{2} 2\varsigma_t(\theta_1+ \nu \theta_2)(I_n -2\varsigma_t(\theta_1+\nu \theta_2))^{-1}2\varsigma_t \theta_2(I_n -2\varsigma_t(\theta_1+\nu \theta_2))^{-1}\right], \\[3pt]g_3(t,\theta_1,\theta_2,\nu) &= \frac{\beta}{2}\text{tr}\left[2\varsigma_t\theta_2(I_n -2\varsigma_t(\theta_1+\nu \theta_2))^{-1} \right],\end{align*}

and $\varsigma_t$ and $\vartheta_t$ as in Proposition 3.1.

We build on the potential approach proposed by Rogers (Reference Rogers1997) whose main idea is to directly specify the state-price density in (4) so that the corresponding instantaneous mortality intensities are positive.Footnote ³ The standard approach proceeds the other way around, that is, to specify the mortality intensities and deduce the state-price density. We briefly summarize the mathematical framework as well as develop our main results on mortality modeling.

First, we define the state-price density as

(16) \begin{align}\zeta_t&\, :\!= \, e^{-\alpha t} (1+ \text{tr}[u_0 v_t])\,, \end{align}

with $\alpha \in \mathbb{R}_+$ , $u_0 = u_1 + u_2$ with $u_1 \in \mathbb{S}_n^{+}$ (the set of positive semi-definite matrices) and $u_2 \in \mathbb{S}_n^{+}$ (so that $u_0 \in \mathbb{S}_n^{+}$ ) and $(v_t)_{t \geq 0}$ a Wishart process. The state-price density can be rewritten as follows. Define the positive function $f \,:\, \mathbb{S}_n^{++} \to \mathbb{R}^+$ such that $f(v) \, :\!= \, 1+\text{tr}[u_0v]$ , then the state-price density writes as $\zeta_t=e^{-\alpha t}f(v_t)$ and is positive as f is. Define $g(v)\, :\!= \,(\alpha - \mathcal{L})f(v)$ and assume that it is a positive function for sufficiently large $\alpha$ . Using Rogers (Reference Rogers1997, Equations (1.2), (1.3), (1.4)), the state-price density allows us to compute $\mathrm{SB}_{0}(t,T)$ given by (4), or any other mortality-related product, if we assume the following change of probability measure

(17) \begin{align}\zeta_t = e^{-\int_0^t (\mu_x(s) + \mu_y(s))ds }\left.\frac{d\mathbb{Q}}{d\mathbb{P}}\right\vert_{\mathcal{F}_t}. \end{align}

Indeed, suppose that one needs to compute the value at time t of $\Pi_t$ that is equal to 1 at time T conditional on the annuitants x and y being alive at that time, then its value is given by

\begin{align*}\Pi_t &= \mathbb{E}_t^{\mathbb{Q}}\left[e^{-\int_t^T (\mu_x(s) + \mu_y(s))ds } \right]= \mathbb{E}_t\left[\frac{\zeta_T}{\zeta_t}\right]\!,\end{align*}

which is known explicitly for a suitable choice of the state-price density $(\zeta_t)_{t\geq 0}$ .Footnote ⁴ The following section shows how the expectations presented in Section 2 can be handled.Footnote ⁵

Our choice of the positive function $ f(v) = 1 + \text{tr}[u_0v] $ is a deliberate aspect of the model specification. The potential approach proposed by Rogers (Reference Rogers1997) necessitates f to be a positive function of the state variable. We selected this particular form of f to ensure that the joint-survival annuity price remains a linear-rational function of the state variable $ v_t $ . It is important to note that alternative choices for f could impact the results discussed later. However, we leave the exploration of these alternatives for future research.

Remark 3.6. The results in this work rely upon the linear-rational structure to obtain the price of the joint-survival annuity and the guaranteed joint-survival annuity option. It is straightforward to see the set of admissible $f \,:\, \mathbb{S}_n^{+} \to \mathbb{R}^+$ that retains the linear-rational structure must take the form

\begin{align*} f(v) = \phi + \operatorname{tr}[\theta_1 v] \end{align*}

where $\phi \in \mathbb{R}$ and $\theta_1$ is a positive semi-definite matrix. Since $\phi$ shifts the level of f(v), it could be difficult to disentangle with the parameters of the Wishart process that control its level. To avoid identification issues, we set $\phi = 1$ .

4. Joint-survival annuity valuation

To illustrate the methodology on the linear-rational Wishart mortality model, let us start with the joint-survival bond that can be explicitly computed as the following proposition shows.

Proposition 4.1. Let $u_1 \in \mathbb{S}_n^+$ and $u_2 \in \mathbb{S}_n^+$ and define $u_0 = u_1 + u_2$ . Assume further that (16) holds. Then on the set $\lbrace \tau_x\gt t, \tau_y\gt t \rbrace $ , the joint-survival bond $\mathrm{SB}(t,T)$ given by (4) is known explicitly as we have

(18) \begin{align}\mathrm{SB}(t,T)&=P(t,T) e^{-\alpha(T-t)}\frac{1 + b_0(T-t) +\text{tr}[a_0(T-t) v_{t}]}{1+\text{tr}[u_0 v_{t}]}, \end{align}

where $ b_0(t)\, :\!= \,b_1(t)+b_2(t)$ and $a_0(t)\, :\!= \,a_1(t)+a_2(t)$ with $(b_1(t),a_1(t))$ and $(b_2(t),a_2(t))$ two sets of functions obtained using Lemma 3.3 for $u_1$ and $u_2$ , respectively.

We should note that the above proposition is identical to the result derived in Da Fonseca (Reference Da Fonseca2024). This is to be expected given the result is independent of the parameterization and the modeling choices we make on the Wishart process.Footnote ⁶ From now on, we will avoid mentioning that the expectations conditional on a given time t are only valid on the sets $\lbrace \tau_x\gt t, \tau_y\gt t \rbrace $ or $\lbrace \tau_x\gt t\rbrace $ , depending on the expectation considered.

According to Rogers (Reference Rogers1997, Equation (1.5)) or Björk (Reference Björk2009, Equation (28.42)) and Cairns (Reference Cairns2004, p. 134), we have

\begin{align*}\mu_x(s)+\mu_y(s) &= \frac{(\alpha - \mathcal{L})f}{f}= \frac{\alpha +\alpha \text{tr}[u_0v_s]-\text{tr}[u_0\omega] - 2 \text{tr}[u_0mv_s]}{1+\text{tr}[u_0v_s]}.\end{align*}

Through our judicious choice of f, we are able to isolate $\alpha$ . Below we define the mortality intensities to be

(19) \begin{align}\mu_x(s) &= \frac{\alpha/2 +\alpha \text{tr}[u_1v_s]-\text{tr}[u_1\omega] - 2 \text{tr}[u_1mv_s]}{1+\text{tr}[u_0v_s]}, \\[-10pt] \nonumber \end{align}

(20) \begin{align}\mu_y(s) &= \frac{\alpha/2 +\alpha \text{tr}[u_2v_s]-\text{tr}[u_2\omega] - 2 \text{tr}[u_2mv_s]}{1+\text{tr}[u_0v_s]}, \\[10pt] \nonumber \end{align}

where we have chosen to apportion the term $\alpha$ equally to both annuitants.Footnote ⁷ From now on, we suppose that there exists $\alpha$ such that $\mu_x(s)$ given by (19) and $\mu_y(s)$ given by (20) are positive. The existence of an $\alpha$ such that the intensities are positive is dependent on the choice of the matrices $u_1, u_2$ and the parameters of the Wishart process. We provide a sufficient requirement in the case when $n=2$ subsequently.

Remark 4.3. Note, we can easily extend our model to the case of an k-joint annuity as follows. First, we should have a dimension of the Wishart process to be at least $n \geq k$ . We then set $u_0 = \sum_{i=1}^k u_i$ , and similar to the above intensities in (19) and (20), we could define the i’th annuitant’s mortality intensity as

\begin{align*} \mu_i(s) = \frac{\alpha/k +\alpha \text{tr}[u_iv_s]-\text{tr}[u_i\omega] - 2 \text{tr}[u_imv_s]}{1+\text{tr}[u_0v_s]}. \end{align*}

In order to understand how the positivity enters into the problem, let us consider the case of $n=2$ (with noting that n could possibly be larger and thus be a multi-factor mortality model). By assumption $u_1\in \mathbb{S}_n^{+}$ and $u_2\in \mathbb{S}_n^{+}$ so that $\text{tr}[u_1\omega]$ , $\text{tr}[u_2\omega]$ , $\text{tr}[u_1v_s]$ , and $\text{tr}[u_2v_s]$ are positive as $\omega \in \mathbb{S}_n^{++}$ while $v_s$ belongs to $\mathbb{S}_n^{++}$ . If $\text{tr}[u_1mv_s]\lt 0$ and $\text{tr}[u_2mv_s]\lt 0$ then as long as $\alpha/2\gt \max\!(\text{tr}[u_1\omega],\text{tr}[u_2\omega])$ , $\mu_x(s)$ and $\mu_y(s)$ are positive. To gain a little bit more of understanding of the positivity problem, it is worth considering the simple case $n=2$ , m diagonal, $u_1=e_{11}$ and $u_2=e_{22}$ then in that particular case $\mu_x(s) $ and $\mu_y(s) $ rewrite as

(21) \begin{align}\mu_x(s) &= \frac{\alpha/2 +\alpha v_{11,s}-\omega_{11} - 2 m_{11}v_{11,s}}{1 + v_{11,s} + v_{22,s}}, \\[-9pt] \nonumber \end{align}

(22) \begin{align}\mu_y(s) &= \frac{\alpha/2 +\alpha v_{22,s}-\omega_{22} - 2 m_{22}v_{22,s}}{1 + v_{11,s} + v_{22,s}}, \\[10pt] \nonumber \end{align}

where the condition $\alpha/2\gt \max\!(\omega_{11},\omega_{22})$ clearly shows that it is sufficient to ensure the positivity of both $\mu_x(s)$ and $\mu_y(s)$ as the matrix m has negative eigenvalues (which means the diagonal terms are negative since m is diagonal). Furthermore, the correlation between these two variables is driven by $d\langle v_{11,.},v_{22,.}\rangle_s$ given by (14), which means it is controlled by the coefficient $(\sigma^2)_{12}$ that can take any sign and is zero if $\sigma_{12}=0$ . As such, $\sigma_{12}$ is the coefficient through which we control the dependency between the two intensities. Analyzing the scaled mortality intensities

(23) \begin{align}{\mu_x^*(s)=(1+\text{tr}[u_0v_s])\mu_x(s), \quad \mu_y^*(s)\, :\!= \,(1+\text{tr}[u_0v_s])\mu_y(s), }\end{align}

we can find their instantaneous quadratic covariation which is given by $d\langle \text{tr}[h_1v_.],\text{tr}[h_2v_.]\rangle_s$ with $h_1=\alpha u_1-2u_1m$ and $h_2=\alpha u_2-2u_2m$ . Using Proposition 3.4, we get

(24) \begin{align} d\langle \mu_x^*(.),\mu_y^*(.)\rangle_s &= d\langle \text{tr}[h_1v_.],\text{tr}[h_2v_.]\rangle_s = \text{tr}[(h_1+h_1^\top) v_s(h_2+h_2^\top)\sigma^2] ds, \end{align}

which demonstrates the impact of the parameters on the covariance between the mortality intensities. Note that when $u_1 = e_{11}$ and $u_2 = e_{22}$ , then the covariation given by (24) will be zero if $\sigma_{12}=0$ .Footnote ⁸ Finally, analyzing the mortality intensities in (21) and (22), given that $\alpha \gt 0$ , $\omega_{ii}\gt 0$ , and $m_{ii} \lt 0$ , we find ${\frac{\partial \mu_x(s)}{\partial v_{11}} \gt 0}$ and $\frac{\partial \mu_x(s)}{\partial v_{22}} \gt 0$ . This provides a natural link between the Wishart states and mortality: as $v_{11,s}$ and $v_{22,s}$ increase, the mortality intensities of x and y also increase, respectively. Thus, the Wishart process can replicate an increasing mortality trend over time in expectation, while remaining stationary.

The expression in (19), combined with the analytical results of Propositions 3.1 as well as the results in Gurland (Reference Gurland1948), enables the computation of the mortality density as the following proposition shows.

Proposition 4.4. Let $\mu_x(T)$ be the mortality intensity of the annuitant x that is given by (19) and denote $c_1=\alpha/2 -\text{tr}[u_1\omega]$ , $h_1 = \alpha u_1-2u_1m$ , then the cumulative distribution function (conditional to time t) of $\mu_x(T)$ is

(25) \begin{align}\mathbb{P}(\mu_x(T) \leq z) = \frac{1}{2}-\frac{1}{\pi}\int_0^{+\infty}\Im \left(\frac{e^{\mathrm{i} s(c_1-z)}\Phi(T-t,\mathrm{i} s(h_1-zu_0),v_t)}{s}\right) ds,\end{align}

where $\Im(.)$ stands for the imaginary part of a complex number and $\Phi(.,.,.)$ is given in Proposition 3.1.

From the above proposition, we can derive the mortality intensity density by simply taking the derivative of the cumulative distribution, which confirms the usefulness of Proposition 3.5, as the following result shows.

Corollary 4.5. Let $\mu_x(T)$ be the mortality intensity of the annuitant x that is given by (19), and thus its cumulative distribution function is given by Proposition 4.4, then the density (conditional to time t) of $\mu_x(T)$ is given by

(26) \begin{align}\frac{d}{dz} \mathbb{P}(\mu_x(T) \leq z) &= -\frac{1}{\pi}\int_0^{+\infty}\Im \left(\frac{(-\mathrm{i}s + g_1+g_2+g_3)e^{\mathrm{i}s(c_1-z)}\Phi(T-t,\mathrm{i}s(h_1-zu_0),v_t)}{s} \right)ds,\end{align}

with $g_j=g_j(T-t,\mathrm{i} sh_1,-\mathrm{i} su_0,z)$ for $j\in \lbrace 1,2,3 \rbrace$ as in Proposition 3.5.

Having an explicit form of the mortality density is of significant practical importance. For example, one can compute the conditional means or higher-order moments of the mortality intensities, which are modeled under the historical measure $\mathbb{P}$ , that can guide in the operational use of the model.

Thanks to (18) the joint-survival bond $\mathrm{SB}(t,T)$ (i.e., (4)) is known and therefore the joint-survival annuity (5) is also known as the next proposition shows.

Proposition 4.7. Let $T_1 \lt \cdots \lt T_N$ be a set of yearly spaced dates such that $t\lt T\lt T_1$ and $T_N=\min\!(x^*-(x+T+ t)-1,y^*-(y+T+t)-1)$ with $x^*$ , resp. $y^*$ , the maximum age the x, resp. y, insured can reach. The joint-survival annuity pays 1 conditional on the survival of both annuitants at those dates, its value at time T is given by

(27) \begin{align}\sum_{i=1}^N \mathrm{SB}(T,T_i)=\frac{b_3(T,T_N)+\text{tr}[a_3(T,T_N)v_T]}{1 + \text{tr}[u_0v_T]}, \end{align}

with

(28) \begin{align}b_3(T,T_N) &\, :\!= \, \sum_{i=1}^N P(T,T_i)e^{-\alpha(T_i-T)} \left( 1 + b_0(T_i-T)\right), \\[-10pt] \nonumber \end{align}

(29) \begin{align}a_3(T,T_N) &\, :\!= \, \sum_{i=1}^N P(T,T_i)e^{-\alpha(T_i-T)} a_0(T_i-T), \\[10pt] \nonumber \end{align}

with $b_0(.), a_0(.)$ given in Proposition 4.1.

Note that the joint-survival annuity (27) also has a linear-rational form, so its density can be computed explicitly (up to a one-dimensional numerical integration).

Proposition 4.8. Let $\sum_{i=1}^N\mathrm{SB}(T,T_i)$ be a joint-survival annuity defined in Proposition 4.7. Then the cumulative distribution function (conditional to time t) of $\sum_{i=1}^N\mathrm{SB}(T,T_i)$ is given by

(30) \begin{align} \mathbb{P}\left(\sum_{i=1}^N\mathrm{SB}(T,T_i) \leq z \right) = \frac{1}{2}-\frac{1}{\pi}\int_0^{+\infty}\Im \left(\frac{e^{\mathrm{i}s(b_3 - z)} \Phi(T-t, \mathrm{i}s(a_3 - z u_0), v_t)}{s}\right) ds,\end{align}

where $a_3$ and $b_3$ are the expressions in Proposition 4.7 with their dependencies dropped on $(T,T_N)$ and $\Phi(.,.,.)$ is given in Proposition 3.1.

Since we have the density of the future annuity value, we can compute certain risk management quantities such as the lower tail that gives the value at risk or the expected shortfall. This sharply contrasts with the mortality model based on the exponential-affine framework for which these quantities are very hard to obtain. Further along this line, an option on the joint-survival annuity, also named the guaranteed joint-survival annuity option, can be priced very efficiently in the linear-rational Wishart mortality model as the next section shows.

5. Guaranteed joint-survival annuity option valuation

With the joint-survival annuity being priced, the next step is to calculate the value of an option on that contract and also to develop option price approximations that simplify the model implementation as well as the risk management of the product. First, we show that an option on the joint-survival annuity is surprisingly simple to price in the linear-rational Wishart mortality model.

Proposition 5.1. The value at time t of a guaranteed joint-survival annuity option (7) with maturity T is given by

(31) \begin{align}\,\bar{\!C}(t,T,T_N) \, :\!= \, P(t,T)e^{-\alpha (T-t)}\frac{\mathbb{E}_t\left[(b_4(T,T_N) + \text{tr}[a_4(T,T_N)v_T])_+\right] }{1+\text{tr}[u_0v_0]},\end{align}

with $b_4(T,T_N)$ a constant and $a_4(T,T_N)$ a matrix such that

(32) \begin{align}b_4(T,T_N) \, :\!= \, b_3(T,T_N) - \frac{1}{g}, \quad a_4(T,T_N) \, :\!= \, a_3(T,T_N) - \frac{1}{g}u_0, \end{align}

where $b_3(.)$ and $a_3(.)$ are given in Proposition 4.7.

As formula (31) shows, the pricing of a guaranteed joint-survival annuity option can be computed using an integration. Indeed, if $Y_T = b_4(T,T_N) + \text{tr}[a_4(T,T_N)v_T] $ , then its characteristic function is given by

(33) \begin{align}\Phi_Y(z) &\, :\!= \,\mathbb{E}_t\left[e^{\mathrm{i} z Y_T} \right]= e^{\mathrm{i} z b_4(T,T_N)}\Phi(T-t,\mathrm{i} z a_4(T,T_N),v_0), \end{align}

with $z\in \mathbb{C}$ , $\mathrm{i}=\sqrt{-1}$ and the function $\Phi(.,.,.)$ given by (11). The expectation in (31) can be computed as

(34) \begin{align}\mathbb{E}_t\left[(Y_T)_+\right]=\frac{1}{\pi}\int_0^{+\infty} \Re\left( \frac{\Phi_Y(z+\mathrm{i}z_i)}{(\mathrm{i}(z+\mathrm{i}z_i))^2} \right) dz\,, \end{align}

where $z_i\lt 0$ and $\Re(.)$ is the real part.

As the result above shows, the valuation of the guaranteed joint-survival annuity option only requires one-dimensional integration. Notice that the result holds whether $\omega =\beta \sigma^2$ (i.e., Bru case) is satisfied or not (see Da Fonseca, Reference Da Fonseca2024 for an implementation in the non-Bru case). In that latter case, computing the moment-generating function involves further numerical integration. Even if the Bru case condition can be a bit too restrictive in practice, it offers nonetheless certain advantages in terms of option price approximations that provide quick alternatives to the exact formula. In what follows, we derive three different approximations to the guaranteed joint-survival annuity option price. These analytic approximations are new and would not be available in Da Fonseca (Reference Da Fonseca2024).

To compute an approximation to the expectation (31), one needs the moments of the variable $Y_T = b_4(T,T_N) + \text{tr}[a_4(T,T_N)v_T]$ , which by definition can be obtained by taking successive derivatives of the moment-generating function. These derivatives can be tedious to perform in the general case, see for example Letac and Massam (Reference Letac and Massam2008). Fortunately, for our Bru case, it is possible to obtain the moments by a direct series expansion.

Proposition 5.2. Suppose that $\omega$ in (8) is such that $\omega=\beta \sigma^2$ with $\beta \in \mathbb{R}$ and $\beta \geq n+1$ and let $Y_T = b_4(T,T_N) + \text{tr}[a_4(T,T_N)v_T] $ with $ b_4(T,T_N)$ and $ a_4(T,T_N)$ (written $ b_4$ and $ a_4$ for simplicity) as in (31). The cumulants of $Y_T$ denoted $(\kappa_j)_{j \in \mathbb{N}_+}$ are given by

(35) \begin{align}\kappa_j &= b_4 \delta_{j1} +\beta (\,j-1)!2^{j-1}\text{tr}[(\varsigma_Ta_4)^{j}]+ j! 2^{j-1} \text{tr}[\vartheta_T^\top(\varsigma_Ta_4)^{\,j}],\end{align}

with $\delta_{j1}=1$ if $j=1$ and 0 otherwise, and $\varsigma_T$ and $\vartheta_T$ given in Proposition 3.1. The moments of $Y_T$ , denoted $(\mu_j)_{j\in \mathbb{N}_+}$ , can be obtained from the cumulants using the relation:

(36) \begin{align}\mu_{j}=\sum _{k=1}^{j}B_{j,k}(\kappa_{1},\ldots ,\kappa_{j-k+1}),\end{align}

where $B_{j,k}$ is the incomplete Bell polynomials.

The cumulants/moments enable the use of Jarrow and Rudd (Reference Jarrow and Rudd1982)’s option price approximation that was also applied in Collin-Dufresne and Goldstein (Reference Collin-Dufresne and Goldstein2002) to the problem of swaption pricing, which is similar to the pricing of an option on a coupon bearing bond or, equivalently, to an option on an annuity. We follow Collin-Dufresne and Goldstein (Reference Collin-Dufresne and Goldstein2002) in our presentation below.

Proposition 5.3. Let $Y_T = b_4(T,T_N) + \text{tr}[a_4(T,T_N)v_T] $ with $ b_4(T,T_N)$ and $ a_4(T,T_N)$ as in (31), and denote the first three cumulants of $Y_T$ , given Proposition 5.2, to be $(\kappa_1,\kappa_2,\kappa_3)$ . Assume that the density of $Y_T$ can be approximated by

(37) \begin{align} \frac{1}{\sqrt{2\pi \kappa_2}} e^{-\frac{(z-\kappa_1)^2}{2\kappa_2}}\left(\sum_{j=0}^3\eta_j (z-\kappa_1)^j\right)\,, \end{align}

with $\lbrace \eta_j;\,j\in 0,\ldots, 3\rbrace$ some constants related to the first three cumulants of $Y_{T}$ . Then

(38) \begin{align}\mathbb{E}_t\left[ \left( Y_{T}\right)_{+}\right] \sim \sum_{j = 0}^3 \eta_j \xi_{j+1}+ \kappa_1\sum_{j = 0 }^{3} \eta_j \xi_{j} \,,\end{align}

where $\lbrace \xi_j;\,j= 0,\ldots, 4 \rbrace$ are some constants that can be computed explicitly and depend on the first three cumulants of $Y_{T}$ .

The above proposition approximates the density of $Y_T$ , the variable involved in the guaranteed joint-survival annuity option, by a perturbation of the Gaussian distribution. Suppose that the eigenvalues of $a_4(T,T_N) + a_4(T,T_N)^\top$ are non-negative, $a_4(T,T_N) + a_4(T,T_N)^\top \in \mathbb{S}_n^{+}$ and hence $\mathrm{tr}[a_4(T,T_N)v_T]\gt 0$ . In this case, the support of the density of $Y_T$ is on the half line and using a Gaussian approximation might not be the most natural choice.

In the Bru case, the Wishart process has a marginal distribution that follows a non-central matrix chi-squared distribution or non-central Wishart distribution. Further can be said when that distribution is projected on a rank-one matrix, it appears in Letac and Massam (Reference Letac and Massam2008, Proposition 3.1) and reads as follows.

The previous proposition gives us an approximation to the characteristic function of the Wishart process for a given T when it is projected onto a rank-one matrix. In fact, in order to apply this proposition, one needs to specify a vector along which the state variable (i.e., the Wishart process) is projected. The following proposition utilizes Proposition 5.4 to obtain an approximation by projecting the state variable in the direction of the eigenvectors associated with the spectral decomposition of $\left(a_4(T,T_N) + a_4(T,T_N)^{\top}\right)/2$ and then summing the distributions obtained from the previous proposition under the assumption of independence. This assumption is inconsistent with the state variable distribution but can still be used to obtain a numerical approximation as detailed below.

Proposition 5.5. Let $Y_T = b_4(T,T_N)+\text{tr}[a_4(T,T_N)v_T]$ of Proposition 5.1 involved in the pricing of the guaranteed annuity option through (33) and (34). Denote the spectral decomposition of

(39) \begin{align}\frac{a_4+a_4^\top }{2} = \sum_{i=1}^n \lambda_i \gamma_i\gamma_i^\top, \end{align}

with $(\gamma_i)_{i=1,\ldots, n}$ some orthonormal vectors and $(\lambda_i)_{i=1,\ldots, n}$ the corresponding (real) eigenvalues (we drop the dependency of $a_4$ and $b_4$ on $(T,T_N)$ ). Rewrite $Y_T$ as

(40) \begin{align}Y_T = b_4+\sum_{i=1}^n \lambda_i \gamma_i^\top v_T \gamma_i, \\[-20pt] \nonumber \end{align}

then according to Proposition 5.4 for each $i=1,\ldots, n$ $\gamma_i^\top v_T \gamma_i/(\varsigma_T)_{ii}$ has a non-central chi-squared distribution with degrees of freedom $\beta$ and non-centrality parameter $(\vartheta_T^\top)_{ii}$ . Then we approximate $Y_T$ with

(41) \begin{align}\tilde{Y}_T = b_4+\sum_{i=1}^n \lambda_i (\varsigma_T)_{ii} \chi_i^2, \\[-20pt] \nonumber \end{align}

where $\left(\chi_i^2(\beta,(\vartheta_T^\top)_{ii})\right)_{i=1,\ldots, n}$ is a set of independent non-central chi-squared random variables (where $\beta$ is the degrees of freedom and each variable has a specific non-centrality parameter $(\vartheta_T^\top)_{ii}$ ). In this case, $\mathbb{E}_t[(Y_T)_+]\sim \mathbb{E}_t[(\,\tilde{\!Y}_T)_+]$ and this latter expectation can be computed using (34) with the characteristic function of the variable $\,\tilde{\!Y}_T $ that depends (linearly) on a generalized chi-squared distribution.

Note that if one eigenvalue in (39) dominates all the others, then the above proposition can be simplified, since the computation of the guaranteed joint-survival annuity option will then only require the integration of the non-central chi-squared distribution, instead of requiring a Fourier transform.

Corollary 5.6. Let $Y_T = b_4(T,T_N)+\text{tr}[a_4(T,T_N)v_T]$ of Proposition 5.1 involved in the pricing of the guaranteed annuity option through (33) and (34). Assume the spectral decomposition (39) and suppose that we have $\lambda_1\geq \lambda_2 \geq \cdots \geq 0 \geq \cdots \geq \lambda_n$ and $\lambda_1\gg |\lambda_j|,\, j = 2,\ldots n$ then approximate $Y_T$ with

(42) \begin{align}\,\tilde{\!Y}_T = b_4+\lambda_1 \gamma_1^\top v_T \gamma_1\,, \end{align}

and according to Proposition 5.4 the scalar variable $\gamma_1^\top v_T \gamma_1/(\varsigma_T)_{11}$ follows a non-central chi-squared distribution with degrees of freedom $\beta$ and non-centrality parameter $(\vartheta_T^\top)_{11}$ .

If $\lambda_1\geq \lambda_2 \geq \cdots \geq 0 \geq \cdots \geq \lambda_n$ and $|\lambda_n| \gg |\lambda_j|,\, j = 1,\ldots n-1$ then approximate $Y_T$ with

(43) \begin{align}\,\tilde{\!Y}_T = b_4 + \lambda_n \gamma_n^\top v_T \gamma_n\,, \end{align}

and according to Proposition 5.4 the scalar variable $\gamma_n^\top v_T \gamma_n/(\varsigma_T)_{11}$ follows a non-central chi-squared distribution distribution with degrees of freedom $\beta$ and non-centrality parameter $(\vartheta_T^\top)_{11}$ .

In Equation (37), we approximate the density of $Y_T$ using a perturbation of the Gaussian distribution, but the marginal distribution of the Wishart process suggests that it is more natural to approximate the density with a perturbation of the (scalar) gamma distribution. Indeed, if $a_4(T,T_N)$ in Proposition 5.2 has positive eigenvalues, then necessarily $b_4\lt 0$ and the support of the density of $Y_T$ is in $[k,\; +\infty)$ with $k\lt 0$ , and if $a_4(T,T_N)$ has negative eigenvalues, then necessarily $b_4\gt 0$ and the support of the density of $Y_T$ is in $(-\infty,\; k]$ for $k\gt 0$ . As such, the Gaussian distribution, whose support is in $\mathbb{R}$ , might not be the most convenient distribution to approximate the density of $Y_T$ . Further to this, taking the limit as $t \to \infty$ in (11) shows that the asymptotic distribution of the Wishart process has a moment-generating function of the form $\frac{1}{\det(I_n - 2\varsigma_t\theta_1)^{\beta/2}}$ , which corresponds to a matrix gamma distribution. Filipović et al. (Reference Filipović, Mayerhofer and Schneider2013) show how to approximate the density of a random variable using a perturbation of the (scalar) gamma distribution, instead of the Gaussian distribution in (37), if the cumulants of the variable are available. As such, Proposition 5.2 allows us to use their results as the following proposition shows.

Proposition 5.7. Let $Y_T = b_4(T,T_N) + \text{tr}[a_4(T,T_N)v_T] $ with $ b_4(T,T_N)$ and $ a_4(T,T_N)$ as in (31) and denote the first three cumulants of $Y_T$ , given Proposition 5.2, to be $(\kappa_1,\kappa_2,\kappa_3)$ . If the eigenvalues of $ a_4(T,T_N)$ are all positive, then we can obtain the following approximation

(44) \begin{align}\mathbb{E}_t\left[ \left( Y_{T}\right)_{+}\right]&= \sum_{j=0}^{3}\sum_{i=0}^{j} \frac{c_j\gamma_{j,i}}{\bar{\beta}} ( \bar{\alpha}+1)_{i+1} Q(\bar{\alpha}+i+2,-\bar{\beta} k) \nonumber \\[3pt]&\quad+ \sum_{j=0}^{3}\sum_{i=0}^{j} c_j\gamma_{j,i}k ( \bar{\alpha}+1)_{i} Q(\bar{\alpha}+i+1,-\bar{\beta} k)\,, \end{align}

with $\lbrace c_j;\, j =0,\ldots, 3 \rbrace$ , $\lbrace \gamma_{j,i};\, j,i =0,\ldots, 3 \rbrace$ , $\bar{\alpha}$ and $\bar{\beta}$ being some parameters that depend on k and the first three moments of $Y_{T}$ , while Q(s,x) stands for the upper regularized incomplete gamma function and $(x)_n$ is the Pochhammer function.

If the eigenvalues of $ a_4(T,T_N)$ are all negative, then we can obtain the following approximation

(45) \begin{align}\mathbb{E}_t\left[ \left( Y_{T}\right)_{+}\right]&= -\sum_{j=0}^{3}\sum_{i=0}^{j} \frac{c_j\gamma_{j,i}}{\bar{\beta}} ( \bar{\alpha}+1)_{i+1} P(\bar{\alpha}+i+2,-\bar{\beta} k) \nonumber\\[3pt]&\quad- \sum_{j=0}^{3}\sum_{i=0}^{j} c_j\gamma_{j,i}k ( \bar{\alpha}+1)_{i} P(\bar{\alpha}+i+1,-\bar{\beta} k)\,, \end{align}

with $\lbrace c_j;\, j =0,\ldots, 3 \rbrace$ , $\lbrace \gamma_{j,i};\, j,i =0,\ldots, 3 \rbrace$ , $\bar{\alpha}$ and $\bar{\beta}$ being some parameters that depend on k and the first three moments of $Y_{T}$ , while P(s, x) stands for the lower regularized incomplete gamma function and $(x)_n$ is the Pochhammer function.

6. Empirical implementation

This section implements the model on the Canadian dependent-lives annuity portfolio canlifins from CASdatasets, the same dataset analyzed by Frees et al. (Reference Frees, Carriere and Valdez1996). We fix the Gompertz marginals at the Frees et al. (Reference Frees, Carriere and Valdez1996) estimates and estimate only the Wishart dependence layer from first-death data. Throughout, we set the risk-free rate to zero ( $r=0$ ) so that in (4) $\mathrm{SB}(0,t)=\mathrm{SB}_0(0,t)$ ; we write SB henceforth, and we write $\mathrm{SB}(0,t\mid v_0)$ when we wish to emphasize the dependence on the contract-specific initial state. Sensitivity analyses, results under alternative model specifications, and a comparison of the three approximation methods (Gaussian, spectral decomposition, and gamma) are reported in the Online Supplementary Appendix. A full empirical calibration that jointly estimates marginals and dependence with a more general loading structure and possibly higher-dimensional factors is deferred to future work.

We reconstruct the observation window for each contract using the annuity guarantee horizon as the censoring boundary, since a censored observation is followed only until the contract expires, not until the raw exit age recorded in the file. For each contract, we define durM and durF as the time to death for the male and female, respectively, or the contract horizon if the corresponding death is not observed. The first-death duration is then $T_1=\min\!(\text{durM},\text{durF})$ , with event indicator $d_1=1$ when at least one death is observed before the contract horizon. We refer to the Supplementary Appendix for data details.

To avoid confounding marginal misspecification with dependence estimation, we fix the Gompertz marginal parameters at the values reported by Frees et al. (Reference Frees, Carriere and Valdez1996). That study estimates two sets of Gompertz parameters: one under an independence assumption (univariate estimation) and one jointly with a Frank copula (bivariate estimation); the two sets differ because the copula-based estimation accounts for the positive dependence between lifetimes. The independence and Wishart specifications use the univariate marginals, while the Frank copula benchmark uses the bivariate marginals from its own joint estimation. Parameter values and corresponding survival curves are reported in the Supplementary Appendix. The conditional survivals $(\,\bar{\!F}_1(t\mid x),\,\bar{\!F}_2(t\mid y))$ and Gompertz entry hazards $(\mu_1^{\text{Gomp}}(x),\mu_2^{\text{Gomp}}(y)$ ) for a given pair of ages (x,y) derived from this specification are used throughout the empirical section. They are given in the Supplementary Appendix.

For each contract $i=1,\ldots,M$ , the data consist of a first-death duration $T_{1,i}$ and a first-death indicator

\begin{align*}d_{1,i}=\begin{cases}1, & \text{if at least one death is observed before the contract horizon},\\[3pt]0, & \text{otherwise (right-censored)}.\end{cases} \end{align*}

A contract is right-censored when the guarantee horizon expires before either death is observed. Either the first death has not yet occurred by $T_{1,i}$ , or its actual timing is unknown. The likelihood contribution for an observed event ( $d_{1,i}=1$ ) is the first-death density $f_1(T_{1,i})$ , while the contribution for a censored observation ( $d_{1,i}=0$ ) is the joint-survival probability $\mathrm{SB}(0,T_{1,i})$ , reflecting the information that both lives were still alive at the censoring time. The standard right-censored first-death log-likelihood is therefore

(46) \begin{align}\ell(\theta)=\sum_{i=1}^M \Bigl[d_{1,i}\log f_1(T_{1,i};\,\theta)+(1-d_{1,i})\log \mathrm{SB}(0,T_{1,i};\,\theta)\Bigr],\end{align}

where $\theta$ is the parameter vector specified below. The likelihood accounts for left truncation through the conditional marginals and for right censoring through the first-death indicator, following the likelihood construction of Frees et al. (Reference Frees, Carriere and Valdez1996, Equations (7)–(13)). Unlike their full bivariate formulation, which uses both individual death times, we use only first-death information, since the joint-survival bond and the first-death density are available in closed form under the Wishart specification. The log-likelihood being specified, the optimal $\theta^*$ satisfies:

(47) \begin{align}\theta^{*}=\operatorname*{argmax}_{\theta}\ell(\theta).\end{align}

Wishart dependence specification: Conditional on the fixed marginals, we estimate the Wishart dependence layer from the first-death likelihood (46). The parameter vector is

(48) \begin{equation}\theta=(\alpha,\beta,m_{11},m_{22},\sigma_{11},\sigma_{22},\rho_\sigma,\lambda). \end{equation}

The parameters $\alpha,\beta,m,\sigma$ are those of the linear-rational Wishart model; the empirical loading coefficient $\lambda$ and the diffusion correlation $\rho_\sigma$ are defined below. The empirical specification takes $n=2$ , keeps m diagonal with negative entries $m_{11},m_{22}\lt 0$ (mean-reverting drift), and parametrizes the loading matrices through a one-parameter common-loading family controlled by $\lambda\in[0,1)$ :

(49) \begin{align}u_1=\begin{pmatrix}1 &\quad \lambda\\[3pt]\lambda &\quad \lambda^2\end{pmatrix}, \qquad u_2=\begin{pmatrix}\lambda^2 &\quad \lambda\\[3pt]\lambda &\quad 1\end{pmatrix}, \qquad u_0=u_1+u_2=\begin{pmatrix}1+\lambda^2 &\quad 2\lambda\\[3pt]2\lambda &\quad 1+\lambda^2\end{pmatrix}.\end{align}

Since $u_1=q_1q_1^\top$ , $u_2=q_2q_2^\top$ with $q_1=(1,\lambda)^\top$ and $q_2=(\lambda,1)^\top$ , each $u_j$ is positive semi-definite and therefore a rank-one matrix. Their sum $u_0$ has eigenvalues $(1+\lambda)^2$ and $(1-\lambda)^2$ , which are strictly positive for $\lambda\in[0,1)$ ; hence, $u_0\in\mathbb{S}_2^{++}$ and therefore rank two. The diffusion matrix is parametrized by $\sigma_{11},\sigma_{22}\gt 0$ and the correlation coefficient $\rho_\sigma\in(-1,1)$ , with off-diagonal entry $\sigma_{12}=\rho_\sigma\sqrt{\sigma_{11}\sigma_{22}}$ so that $\sigma$ is positive definite.

This parametrization is a parsimonious alternative to a fully unrestricted loading structure: the model allows a common factor to enter both intensities but remains simple enough to preserve closed-form first-death quantities while permitting the analytical recovery of the contract-specific initial state described below.

We restrict the initial state to a diagonal matrix,

(50) \begin{align}v_0=\begin{pmatrix}v_{11,0} &\quad 0\\[3pt]0 &\quad v_{22,0}\end{pmatrix},\end{align}

where the diagonal entries $v_{11,0}$ and $v_{22,0}$ are contract-specific and determined by matching the model’s time-zero mortality intensities to the fixed Gompertz entry hazards for the contract’s entry ages (x,y); the inversion procedure is described hereafter. The off-diagonal component $v_{12,0}$ is set to zero at entry, but the common loading $\lambda$ and the diffusion correlation $\rho_\sigma$ still affect the subsequent dynamics and hence the first-death survival law.

Specifying $m, v_0$ diagonal and rank-one loading matrices (i.e., $u_1,u_2$ ) represents a parsimonious model parametrization that enables identification from first-death data alone. Note, however, that there are three channels through which the model generates dependence between the mortality intensities (19) and (20): (i) structural dependence through the shared denominator $1+\text{tr}[u_0 v_s]$ , which is present even when $\lambda=0$ ; (ii) the common loading controlled by $\lambda$ , which allows each life’s intensity to depend on the other life’s state component through the off-diagonal entries of $u_1$ and $u_2$ ; and (iii) the diffusion correlation controlled by $\sigma_{12}$ , which governs the covariation between the scaled intensities (24) through the instantaneous covariation between $v_{11,t}$ and $v_{22,t}$ given by (14). Any relaxation of these restrictions introduces parameters that the first-death likelihood cannot separately identify. A non-diagonal m would add off-diagonal drift terms that are confounded with the loading structure; a non-diagonal $v_0$ would require knowledge of the initial cross-state $v_{12}$ , which is not pinned down by the two entry hazards; and an unrestricted loading choice $u_1,u_2\in\mathbb{S}_2^+$ would introduce six free entries, but the joint-survival bond depends on $u_1$ and $u_2$ only through $u_0=u_1+u_2$ and the trace combinations $\text{tr}[u_j\omega]$ and $\text{tr}[u_j m v]$ .

Under diagonal m, the matrix $\mathsf{A}$ in Lemma 3.3, used in the pricing formulas of Proposition 4.1, is diagonal, so its exponential is a diagonal matrix with the (scalar) exponential of the diagonal elements of $\mathsf{A}$ on its diagonal. Because $v_0$ is diagonal, the off-diagonal entries of $a_0(t)$ do not contribute to $\text{tr}[a_0(t)v_0]$ and therefore to the numerator of a survival bond (18). Explicit expression of the survival bond, which is reported in Supplementary Appendix, clarifies the model identification through survival bonds. The parameter $\lambda$ appears in the numerator and the denominator of a survival bond, while the diffusion correlation $\rho_\sigma$ affects the sign and magnitude of $\omega_{12}=\beta(\sigma^2)_{12}$ , which is shared by the two survival bonds. When $\lambda=0$ and $\rho_\sigma=0$ , the cross term vanishes and the specification collapses to two independent one-dimensional Wishart (i.e., Cox–Ingersoll–Ross) processes. As explained in Section 4, a sufficient condition for non-negative intensities is $\alpha/2\gt \max\!(\text{tr}[u_1\omega],\text{tr}[u_2\omega])$ this is enforced during estimation.

The model estimation loops over a two-stage optimization procedure. One of the stages consists of a given parameter vector $\theta$ to estimate $v_0$ . To this end, we assume that the global parameter vector $\theta$ is shared across all contracts, and each contract $i=1,\ldots,M$ with entry ages $(x_i,y_i)$ requires its own initial state $v_0^{(i)}=\text{diag}(v_{11,0}^{(i)},v_{22,0}^{(i)})$ so that the model’s time-zero mortality intensities match the Gompertz entry hazards denoted $\mu_1^{\text{Gomp}}(x_i)$ and $\mu_2^{\text{Gomp}}(y_i)$ . It means that setting $t=0$ in (19) and (20) along with the vector parameter $\theta$ , the matching conditions for a generic contract with entry ages $(x_i,y_i)$ read

(51) \begin{align}\mu_x(0) &= \mu_1^{\text{Gomp}}(x_i),\, \quad \mu_y(0) = \mu_2^{\text{Gomp}}(y_i), \end{align}

where the intensities $\mu_x(0),\mu_y(0)$ depend on i through the initial state variable value $v_0^{(i)}$ . Since these intensities are linear-rational in the state variable, the above equalities can be rewritten as

(52) \begin{align}\mathbf{A}_i\begin{pmatrix}v^{(i)}_{11,0}\\[3pt] v^{(i)}_{22,0}\end{pmatrix}= \mathbf{b}_i, \end{align}

with $\mathbf{A}_i$ a matrix and $\mathbf{b}_i$ a vector (we refer to the Supplementary Appendix for the expressions). The system (52) has a unique solution whenever $\det(\mathbf{A}_i)\neq 0$ . Since $\theta$ is global, the coefficient matrix $\mathbf{A}_i$ changes from contract to contract only through the Gompertz entry hazards $\mu_1^{\text{Gomp}}(x_i)$ and $\mu_2^{\text{Gomp}}(y_i)$ ; the solution for each contract is computed by Cramer’s rule. The parameter vector is accepted only if $v^{(i)}_{11,0}\gt 0$ and $v^{(i)}_{22,0}\gt 0$ for every contract in the sample; this positivity requirement is enforced during estimation. This stage is combined with a second stage that consists of a set of fixed values $ (v_0^{(i)})_{i=1,\ldots,M}$ to solve (47). As a result, it remains to clarify how $\ell(\theta)$ of (46) is determined in the Wishart model.

The survival bond $\mathrm{SB}(0,t\mid v_0)$ of (18) is explicitly known and equals the joint-survival probability so that the first-death density is obtained by differentiation:

(53) \begin{align}f_1(t\mid v_0)=-\frac{d}{dt}\mathrm{SB}(0,t\mid v_0)=\mathrm{SB}(0,t|v_0)(\mu_x(0)+\mu_y(0)),\end{align}

and is positive whenever the intensities (19) and (20) are non-negative which is guaranteed by the choice of the model parameters. As a result, both the joint-survival bond and the first-death density (53) are available in closed form, so no numerical integration or simulation is required inside the likelihood evaluation (46) where each component of the sum is contract-level dependent through its initial state $v_0^{(i)}$ (which is related to the pair $(x_i,y_i)$ through $(\mu_1^{\text{Gomp}}(x_i),\mu_2^{\text{Gomp}}(y_i))$ ) and depends on the global parameter vector $\theta$ through the inversion (52) and the expression of the survival bond.

The optimization is carried out over eight unconstrained variables, each mapped to its natural parameter via a smooth bijection (exponential, hyperbolic tangent, or a shift thereof) that enforces the relevant domain constraint. The negative log-likelihood (NLL) $-\ell(\theta)$ is minimized by a deterministic multi-start procedure. The parameter transforms and the optimization protocol are described in the Supplementary Appendix.

We compare the linear-rational Wishart model to two alternatives, the independent model and the fixed-parameter Frank copula model. The comparison is performed on the same first-death data and fixed Gompertz marginals, all evaluated under the standard first-death log-likelihood (46). Under independence, the joint survival is the product of the Gompertz conditional survivals

(54) \begin{align}\mathrm{SB}^{\text{ind}}(0,t)=\,\bar{\!F}_1(t\mid x)\, \,\bar{\!F}_2(t\mid y), \end{align}

with $(\,\bar{\!F}_1(t\mid x),\,\bar{\!F}_2(t\mid y))$ being the single-life conditional survival functions, while the first-death density is

(55) \begin{align}f_1^{\text{ind}}(t)=\mathrm{SB}^{\text{ind}}(0,t)\bigl(\mu_1^{\text{Gomp}}(x+t)+\mu_2^{\text{Gomp}}(y+t)\bigr). \end{align}

This model has no free dependence parameters. For the Frank copula, we follow Frees et al. (Reference Frees, Carriere and Valdez1996) and use $C_\delta(u,v)=-(1/\delta)\log\!\bigl(1+(e^{-\delta u}-1)(e^{-\delta v}-1)/(e^{-\delta}-1)\bigr)$ with the literature parameter $\delta=3.367$ . The joint survival is

(56) \begin{align}\mathrm{SB}^{\text{Frank}}(0,t)=C_\delta \bigl(\,\bar{\!F}_1(t\mid x),\,\,\bar{\!F}_2(t\mid y)\bigr),\end{align}

and the first-death density is obtained by differentiating $-\mathrm{SB}^{\text{Frank}}$ with respect to t. Because the copula parameter is fixed at its literature value, the Frank benchmark also has no free dependence parameters in this comparison. For the Wishart-based model, the specification is the one mentioned above with the all eight parameters of (48) free.

Table 1.

Estimated Wishart dependence parameters from the first-death likelihood.

Table 1. Long description

A table with two panels, Panel A and Panel B, displaying estimated Wishart dependence parameters and initial state variables. Panel A: Wishart parameters. This panel contains a table with two rows and eight columns. The columns are labeled Parameter and Value. The parameters listed are alpha, beta, m11, m22, sigma11, sigma22, rho_sigma, and lambda. The corresponding values are 0.033, 3.073, -0.370, -0.237, 0.070, 0.072, 0.127, and 0.256. Panel B: Initial state variable. This panel contains a table with six rows and five columns. The columns are labeled Age pair, v11_0, and v22_0. The age pairs listed are (60,60), (65,60), (70,65), (75,70), and (80,75). The corresponding values for v11_0 are 0.0117, 0.0179, 0.0283, 0.0463, and 0.0783. The corresponding values for v22_0 are 0.0090, 0.0084, 0.0115, 0.0177, and 0.0303.

Note: Panel A, estimates for (48) by solving (47), the maximization of the right-censored first-death log-likelihood (46) with fixed Gompertz marginals from Frees et al. (Reference Frees, Carriere and Valdez1996). The shared-loading coefficient $\lambda$ and the diffusion correlation $\rho_\sigma$ govern the dependence between the intensities. Panel B, solutions to the system (52) for different pairs of ages.

Estimation results: Table 1 reports the estimated coefficients, while the NLLs for the three models are Independence $8114.56$ , Frank $8129.35$ , and Wishart $8322.38$ (with further details in the Supplementary Appendix). Among the models compared, independence achieves the lowest NLL, followed by the Frank copula and the Wishart specification. Empirically, the Wishart fit identifies non-trivial off-diagonal dependence: $\widehat{\lambda}=0.256$ and $\widehat{\rho}_\sigma=0.127$ . To quantify the dependence implied by these estimates, we evaluate the instantaneous correlation between the scaled mortality intensities (23) at the stationary mean of the Wishart process. Indeed, thanks to (24), this instantaneous correlation is

(57) \begin{align}\text{Corr}(d\mu_x^*(s), d\mu_y^*(s)) = \frac{\text{tr}\bigl[(h_1+h_1^\top)v_s(h_2+h_2^\top)\sigma^2\bigr]}{\sqrt{\text{tr}\bigl[(h_1+h_1^\top)v_s(h_1+h_1^\top)\sigma^2\bigr]\text{tr}\bigl[(h_2+h_2^\top)v_s(h_2+h_2^\top)\sigma^2\bigr]}},\end{align}

but because the initial state is diagonal with $v_{12}=0$ , evaluating it at $v_0$ yields a correlation that does not reflect the dependence generated by the subsequent Wishart dynamics. Instead, we compute it at the stationary mean $\bar{v}_\infty\, :\!= \,\lim_{t\to\infty}\mathbb{E}[v_t]$ , which solves $\text{vec}(\bar{v}_\infty)=-\mathsf{A}^{-1}\text{vec}(\omega)$ with $\mathsf{A}$ as in Lemma 3.3, that is given by

\begin{align*}\bar{v}_\infty=\begin{pmatrix}{\color{black}0.0209} &\quad {\color{black}0.0065}\\[3pt]{\color{black}0.0065} &\quad {\color{black}0.0338}\end{pmatrix}. \end{align*}

With the estimated parameters reported in Table 1, evaluating (57) at $\bar{v}_\infty$ yields an instantaneous correlation of approximately ${\color{black}0.430}$ , indicating that the common-loading Wishart specification implies a stationary instantaneous correlation of about ${\color{black}43\%}$ between the scaled mortality intensities. This value falls within the interval $(0.41,0.56)$ reported by Frees et al. (Reference Frees, Carriere and Valdez1996, p. 237) for the correlation between annuitants in last-survivor contracts, suggesting that even the parsimonious common-loading family captures the positive dependence documented in the data. Extensions to a non-diagonal initial state allowing for a time-varying correlation profile from inception are left for future work.

The likelihood ranking still requires careful interpretation, as the comparison is structurally asymmetric. The independence and Frank copula benchmarks evaluate the Gompertz conditional survival function $\,\bar{\!F}_j(t\mid x)$ at every horizon t, so they benefit from the exact Gompertz marginal shape at all times. The Wishart specification, by contrast, matches the Gompertz entry hazards only at $t=0$ through the inversion (52); for $t\gt 0$ , the Wishart dynamics govern the mortality intensity, and the implied survival curve necessarily diverges from Gompertz. The first-death likelihood therefore penalizes the Wishart primarily for its marginal shape rather than for any deficiency in its dependence specification, in line with the evidence of De Giovanni et al. (Reference De Giovanni, Pirra and Viviano2025) that, on this same portfolio, marginal flexibility rather than the dependence layer governs goodness-of-fit.Footnote ⁹ Furthermore, since the Wishart model is not nested in the other specifications, comparing likelihood values convey little information. It might be more informative to simply compare the survival bond pricing errors. As the joint-life annuity under the fitted specification is computed by summation of the closed-form survival bonds, for annual payments on a horizon of N years with zero interest rates, the joint-life annuity (i.e., (27)) at entry is

(58) \begin{align}\mathcal{A}_{xy}=\sum_{k=1}^N \mathrm{SB}(0,k\mid v_0), \end{align}

with SB given by (18) and the contract-specific $v_0$ obtained from the entry-hazard inversion (52).

The fitted models imply different annuity values when evaluated under fixed Gompertz marginal inputs from Frees et al. (Reference Frees, Carriere and Valdez1996), with the independence and Wishart specifications using the univariate marginals and the Frank benchmark using the bivariate marginals from the same study. Table 2 reports joint-life annuity values for annual payments, zero interest rates, and a five-year horizon equal to the observed contract support, for the age pairs (60, 60), (65, 60), (70, 65), (75, 70), and (80, 75). Thus, the annuity comparison reflects differences in both the marginal inputs and the joint-survival specification. At entry ages (60, 60), the Wishart joint-life annuity is $4.68$ , compared with $4.83$ under independence (a decrease of $0.15$ , or $\sim$ $3.0\%$ ); at (65, 60) the gap narrows to $-0.10$ ( $4.66$ vs. $4.75$ ). The sign reverses for older entrants: at (75,70) the Wishart value is $4.52$ against $4.33$ under independence (an increase of $0.19$ , or $\sim$ $4.5\%$ ), rising to a $0.44$ difference ( $4.37$ vs. $3.93$ , or $\sim$ 11%) at (80,75). The pattern reflects the structural divergence between the Wishart and Gompertz survival curves: at younger ages, where the Wishart dynamics match the Gompertz marginals closely, the model implies slightly lower joint-survival probabilities; at older ages, where the Gompertz hazard rises exponentially but the mean-reverting Wishart intensity reverts to a finite long-run level, the model implies higher joint survival and hence a larger joint-life annuity.

Table 2.

Joint-life annuity values $\mathcal{A}_{xy}$ across benchmark models.

Table 2. Long description

A table comparing joint-life annuity values across different age pairs and models. The table has five rows and five columns. The columns are labeled with age pairs: (60,60), (65,60), (70,65), (75,70), and (80,75). The rows are labeled with different models: Ind., Frank, Wishart, and Δ. The values in the table are as follows: Row 1: Ind., 4.83, 4.75, 4.59, 4.33, 3.93. Row 2: Frank, 4.84, 4.76, 4.60, 4.36, 4.00. Row 3: Wishart, 4.68, 4.66, 4.61, 4.52, 4.37. Row 4: Δ, 0.01, 0.19, 0.44. The table shows the joint-life annuity values for annual payments, zero interest rates, and a five-year horizon equal to the observed contract support. The values reflect differences in both the marginal inputs and the joint-survival specification. The table highlights how the annuity values change across different age pairs and models, with notable trends and comparisons visible in the data.

Note: Annual payments, zero interest rate ( $r=0$ ), and a 5-year horizon equal to the observed contract support. Independence and Frank values use the product and copula joint survivals, respectively, with fixed Gompertz marginals from Frees et al. (Reference Frees, Carriere and Valdez1996). Wishart values are computed from the closed-form survival bond of Proposition 4.1 via (76), with the contract-specific initial state from (52). $\Delta=\text{Wishart}-\text{Independence}$ .

Remark 6.2. The last-survivor annuity requires model-consistent single-life survivals. Under the Wishart specification, the Radon–Nikodým decomposition (17) yields

\begin{align*}\mathbb{E}_t^{\mathbb{Q}}\left[e^{-\int_t^T \mu_x(s)ds}\right] = \mathbb{E}_t^{\mathbb{P}}\left[\frac{\zeta_T}{\zeta_t}e^{\int_t^T \mu_y(s)ds}\right], \end{align*}

which can be estimated by Monte Carlo simulation under $\mathbb{P}$ ; a symmetric expression holds for y. Computing both model-implied single-life annuities by simulation and combining them with the closed-form joint-life annuity via inclusion–exclusion gives a fully model-consistent last-survivor annuity. This extension is deferred to future work.

Guaranteed annuity option pricing: Although the Wishart specification ranks below the copula benchmarks on the first-death likelihood, that comparison evaluates only the static joint-survival function, the one object common to both frameworks. The guaranteed joint-survival annuity option $\,\bar{\!C}(0,T,T_N)$ (see (6)) illustrates a qualitatively different output that only the linear-rational Wishart model can produce, and which lies entirely beyond the scope of any static copula and the standard exponential-affine framework of Dahl (Reference Dahl2004) and Biffis (Reference Biffis2005) is also numerically unmanageable, see Biffis and Millossovich (Reference Biffis and Millossovich2006). Note that the advantages of our proposed framework go well beyond the problem of pricing guaranteed annuity options. Indeed, not only can this product be priced, but risk measures such as value at risk and expected shortfall, to name just a few, can also be explicitly computed for a joint-survival bond. These measures require a dynamic evolution of joint mortality intensities, which the copula approach does not provide, while the multiple-payment feature of survival annuities makes their treatment problematic within the exponential-affine framework. As a result, even at the level of risk management of the underlying asset of the guaranteed annuity option (GAO), our framework is considerably more operational than existing alternatives.

The guaranteed annuity option pays $(\mathcal{A}_{xy}(T,T_N)-1/g)_+$ at maturity T, where $\mathcal{A}_{xy}(T,T_N)$ is the joint-life annuity (58) (or Proposition 4.7) with payment dates $T_1,\ldots,T_N$ and g is the guarantee rate. The specification adopted here is summarized in the appendix but report here the main characteristics: we set the option maturity equal to the observed contract horizon, $T=5$ years, with annual annuity payments at $T_i=i$ for $i=1,\ldots,5$ , and set the guarantee rate at the at-the-money level $g=1/\mathcal{A}_{xy}$ so that the option is exercised whenever the realized annuity value at maturity exceeds its time-zero expectation. By Proposition 5.1, the Wishart framework prices $\,\bar{\!C}(0,T,T_N)$ in semi-closed form via Fourier inversion, and three closed-form alternatives are available: the Gaussian approximation (Proposition 5.3), the spectral-decomposition approximation (Proposition 5.5), and the gamma approximation (Proposition 5.7).

Table 3 reports the resulting prices and the percent relative error of each approximation against the Fourier benchmark. The at-the-money guarantee rates range from $g=0.2137$ at (60,60) to $g=0.2290$ at (80, 75), reflecting the decrease in the underlying annuity value $\mathcal{A}_{xy}$ (from $4.68$ to $4.37$ ) as entry ages rise and mortality risk increases. The Fourier option price rises sharply with entry age: it is $0.002948$ at (60, 60), $0.006358$ at (65, 60), $0.019582$ at (70, 65), $0.057577$ at (75, 70), and $0.147511$ at (80, 75), a 53-fold increase from the youngest to the oldest pair. The steep gradient has a clear actuarial explanation: at younger entry ages the five-year mortality risk is small and the annuity value at maturity is tightly concentrated around its expectation, so the at-the-money option is deep out of the money and nearly worthless; at older ages, the higher mortality hazard generates substantially more variance in the future annuity value, pushing a larger mass of the distribution above the strike $1/g$ and making the guarantee economically material. For instance, at (80, 75) the option price $0.147511$ represents $3.37\%$ of the underlying annuity value of $4.37$ , a non-negligible component of the annuity liability that an insurer must reserve against. Neither independence nor the Frank copula provides any mechanism to compute this price. All four Wishart-based methods are evaluated from the single estimated parameter vector $\theta$ . The spectral approximation is the most accurate across the age grid, with a relative error below $3.5\%$ for all age pairs with entry ages (65,60) and above and below $0.3\%$ at the oldest pair. The gamma approximation performs much better than the Gaussian approximation. This is to be expected: the stationary distribution of the Wishart process is a matrix gamma distribution, so the scalar gamma expansion of Proposition 5.7 inherits the correct distributional shape and captures the tail behavior of the annuity value more faithfully than the symmetric Gaussian expansion. The Gaussian approximation is comparably accurate at older ages (below 5% for (75, 70) and (80, 75)) but overestimates substantially when the option is deep out of the money, precisely the regime where the Gaussian tail overshoot relative to the gamma tail is most pronounced. The gamma approximation performs well for all the pairs with larger errors at younger ages. Overall, these relative errors are largest at the youngest age pairs, where the option price itself is very small; the absolute pricing differences across all methods remain modest throughout for longer ages.

Table 3.

Guaranteed joint-survival annuity option prices $\,\bar{\!C}(0,T,T_N)$ from the fitted Wishart model.

Table 3. Long description

A table comparing guaranteed joint-survival annuity option prices from the fitted Wishart model across different age pairs and methods. The table has six columns and five rows, including the header row. The columns are labeled Age pair, g, Fourier, Gaussian, % err, Spectral, % err, and Gamma, % err. The row labels are age pairs ranging from (60,60) to (80,75). The values in the table represent the option prices and percent relative errors for each method against the Fourier benchmark. Row 1: Age pair, (60,60); g, 4.6789; Fourier, 0.2137; Gaussian, 0.002948; % err, 0.004198; Spectral, 0.002827; % err, 0.003518; Gamma, % err, 0.003518. Row 2: Age pair, (65,60); g, 4.6563; Fourier, 0.2148; Gaussian, 0.006358; % err, 0.008288; Spectral, 0.006145; % err, 0.007026; Gamma, % err, 0.007026. Row 3: Age pair, (70,65); g, 4.6052; Fourier, 0.2171; Gaussian, 0.019582; % err, 0.021990; Spectral, 0.019171; % err, 0.020074; Gamma, % err, 0.020074. Row 4: Age pair, (75,70); g, 4.5178; Fourier, 0.2213; Gaussian, 0.057577; % err, 0.059918; Spectral, 0.057009; % err, 0.057373; Gamma, % err, 0.057373. Row 5: Age pair, (80,75); g, 4.3677; Fourier, 0.2290; Gaussian, 0.147511; % err, 0.148798; Spectral, 0.147072; % err, 0.147114; Gamma, % err, 0.147114. The table shows that the Fourier option price rises sharply with entry age, reflecting the increase in mortality risk and the underlying annuity value. The spectral approximation is the most accurate across the age grid, with a relative error below 5 percent for all age pairs with entry ages (65,60) and above. The gamma approximation performs better than the Gaussian approximation, capturing the tail behavior of the annuity value more faithfully.

Note: The option payoff is $(\mathcal{A}_{xy}(T,T_N)-1/g)_+$ with parameters from Table 8. Fourier: semi-closed-form Fourier inversion (Proposition 5.1); Gaussian: closed-form normal approximation (Proposition 5.3); Spectral: eigenvalue-projection approximation (Proposition 5.5); Gamma: closed-form gamma approximation (Proposition 5.7). The at-the-money guarantee rate is $g=1/\mathcal{A}_{xy}$ . Percent errors are relative to the Fourier benchmark. All methods use the estimated Wishart parameter vector $\theta$ from Table 1.

In the Supplementary Appendix, a figure plots the guaranteed annuity option price as a function of the guarantee rate g for the entry-age pair (65, 60). It can be checked, at least visually, that the four pricing methods produce qualitatively similar curves over this range, with the spectral approximation tracking the Fourier benchmark most closely. The option price increases with g (since a higher guarantee rate means a more valuable floor) and is always non-negative under the Fourier, spectral, and gamma methods.

This pricing exercise demonstrates the end-to-end tractability that distinguishes the linear-rational Wishart framework from both static copula specifications and exponential-affine intensity models. The same parameter vector $\theta$ that enters the first-death likelihood is fed directly into the guaranteed annuity option price formula, the cumulant expansions, and the closed-form approximations, all without simulation. The first-death likelihood comparison therefore captures only one dimension of model performance; on the pricing and risk-management quantities that motivate the theoretical development, the linear-rational Wishart model is the only specification among those considered that delivers quantitative answers.

This empirical section is best understood as a proof-of-concept implementation that deliberately adopts the simplest admissible specification at every stage. Each of these choices admits natural generalizations. A richer function f, which relates the state-price density and the state variable, could improve the marginal fit of the model-implied survival curves but would require re-deriving the intensity expressions, the positivity conditions, and all downstream pricing formulae. Relaxing the diagonal constraint on $v_0$ introduces an off-diagonal initial correlation but necessitates new sufficient conditions for non-negative intensities that go beyond the analysis developed here. An unrestricted loading structure for $u_1,u_2$ raises identifiability issues that must be resolved before estimation is carried out. Finally, jointly estimating the marginal and dependence parameters, rather than conditioning on a literature benchmark, requires a fundamentally different calibration protocol and careful treatment of the resulting identification constraints. Each of these generalizations involves substantial additional analytical development, and together they constitute a substantial piece of work in its own right, which we leave for future investigation.