2. Risk factors

Our objective is to demonstrate the effectiveness of PINNs in valuing GMABs and establishing investment strategies for them. To validate our methodology, we compare the exact GMAB prices to the estimates generated by the neural network. To facilitate this comparison, we consider Brownian financial and biometric risk factors and derive a closed-form expression for the GMAB price. An alternative validation approach would involve evaluating GMAB prices through MC simulations, but this method is computationally intensive and less accurate. This section introduces the dynamics of risk factors, while Section 3 provides analytical expressions for bond prices and survival probabilities. Finally, we obtain the analytical value of the GMAB in Section 4.

The GMABs provide policyholders with a promise that the value of their investment, conditionally to their survival, will reach a minimum predetermined amount at expiry, regardless of how the underlying investments perform. The value of GMABs depends both on biometrical and financial risk factors that we detail in this section. We consider an hybrid market in which we invest in cash, stocks, interest rate bonds, and mortality bonds. The stock indice and the interest rate are, respectively, denoted by $\left (S_{t}\right )_{t\geq 0}$ , $\left (r_{t}\right )_{t\geq 0}$ . The insured is $x_{\mu }$ year old, and the reference force of mortality is denoted by $\left (\mu _{x+t}\right )_{t\geq 0}$ . The insurer can also invest in mortality bonds written on a reference population of age $x_{\lambda }$ and mortality rates $\left (\lambda _{x_{\lambda }+t}\right )_{t\geq 0}$ .

The financial and biometric processes are defined on a probability space $\left ({\Omega },\mathcal{F},\mathbb{P}\right )$ associated to four independant Brownian motions under $\mathbb{P}$ , ${\tilde{\boldsymbol{W}}}_{t}=\left (\tilde{W}_{t}^{(1)},\tilde{W}_{t}^{(2)},\tilde{W}_{t}^{(3)},\tilde{W}_{t}^{(4)}\right )_{t\geq 0}^{\top }$ . The state variables $\left (S_{t},r_{t},\mu _{x_{\mu }+t},\lambda _{x_{\lambda }+t}\right )$ are driven by the following SDEs:

(1) \begin{eqnarray} d\left (\begin{array}{c} S_{t}\\ r_{t}\\ \mu _{x_{\mu }+t}\\ \lambda _{x_{\lambda }+t} \end{array}\right ) & = & \left (\begin{array}{c} \left (r_{t}+\nu _{S}\sigma _{S}\right )\,S_{t}\\ \kappa _{r}\left (\gamma _{r}(t)-\nu _{r}\frac{\sigma _{r}}{\kappa _{\kappa }}-r_{t}\right )\\ \kappa _{\mu }\left (\gamma _{\mu }(t)-\nu _{\mu }\frac{\sigma _{\mu }(t)}{\kappa _{\mu }}-\mu _{x_{\mu }+t}\right )\\ \kappa _{\lambda }\left (\gamma _{\lambda }(t)-\nu _{\lambda }\frac{\sigma _{\lambda }(t)}{\kappa _{\lambda }}-\lambda _{x_{\lambda }+t}\right ) \end{array}\right )dt+\left (\begin{array}{c} S_{t}\sigma _{S}\Sigma _{S}^{\top }\\ \sigma _{r}\Sigma _{r}^{\top }\\ \sigma _{\mu }(t)\Sigma _{\mu }^{\top }\\ \sigma _{\lambda }(t)\Sigma _{\lambda }^{\top } \end{array}\right )d{\tilde{\boldsymbol{W}}}_{t}\,, \end{eqnarray}

where $\kappa _{r}$ , $\kappa _{\mu }$ , $\kappa _{\lambda }$ , $\sigma _{S}$ , and $\sigma _{r}$ belong to $\mathbb{R}^{+}$ , whereas $\gamma _{r}(t)$ , $\gamma _{\mu }(t)$ , $\gamma _{\lambda }(t)$ , $\sigma _{\mu }(t)$ , and $\sigma _{\lambda }(t)$ are positive functions of time. $\gamma _{r}(t)$ $\gamma _{\mu }(t)$ , and $\gamma _{\lambda }(t)$ are fitted to term structures of interest and mortality rates. Details about the estimation procedure are provided in Appendix. Notice that $\gamma _{\mu }(t)$ , $\gamma _{\lambda }(t)$ , $\sigma _{\mu }(t)$ , and $\sigma _{\lambda }(t)$ depend on the initial age of the insured and of the hedging population; $x_{\mu },x_{\lambda }\in [0,x_{max})$ , $x_{max}\gt 0$ .

In our framework, mortality rates have a mean-reverting Gaussian dynamic, which offers a good trade-off between complexity and analytical tractability. Milevsky & Promislow (Reference Milevsky and Promislow2001) were among the firsts to consider mean-reverting stochastic processes for modeling mortality. Luciano & Vigna (Reference Luciano and Vigna2005) developed a jump-diffusion affine model for modeling rates and showed that constant mean reverting process are not adapted for describing the mortality. Luciano & Vigna (Reference Luciano and Vigna2008) calibrated various time homogeneous mean-reverting models to different generations in the UK population and investigate their empirical appropriateness. Jevtic et al. (Reference Jevtic, Luciano and Vigna2013) studied and calibrated a cohort-based model which captures the characteristics of a mortality surface with a continuous-time factor approach. Zeddouk & Devolder (Reference Zeddouk and Devolder2020) modeled the force of mortality with mean reverting affine processes, where the level of reversion is time dependent. Hainaut (Reference Hainaut2023) proposed a mortality model based on a mean reverting random field indexed by time and age.

In Equation (1), the parameters $\nu _{S}$ , $\nu _{r}$ , $\nu _{\mu }$ , and $\nu _{\lambda }$ tune the risk premiums $\nu _{S}\sigma _{S}$ , $-\nu _{r}\sigma _{r}$ , $-\nu _{\mu }\sigma _{\mu }(t)$ , $-\nu _{\lambda }\sigma _{\lambda }(t)$ of processes. Furthermore, we assume that the standard deviation of mortality rates is related to age through the relations $\sigma _{\mu }(t)=\alpha _{\mu }e^{\beta _{\mu }(x_{\mu }+t)}$ and $\sigma _{\lambda }(t)=\alpha _{\lambda }e^{\beta _{\lambda }(x_{\lambda }+t)}$ . $\Sigma _{S}$ , $\Sigma _{r}$ , $\Sigma _{\mu }$ , and $\Sigma _{\lambda }$ are vectors such that $\Sigma =\left (\Sigma _{S}^{\top },\Sigma _{r}^{\top }\Sigma _{\mu }^{\top },\Sigma _{\lambda }^{\top }\right )$ is the (upper) Choleski decomposition of the correlation matrix:

\begin{eqnarray*} \Sigma & =\left (\begin{array}{c} \Sigma _{S}^{\top }\\[5pt] \Sigma _{r}^{\top }\\[5pt] \Sigma _{\mu }^{\top }\\[5pt] \Sigma _{\lambda }^{\top } \end{array}\right )= & \left (\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} \epsilon _{SS} & \epsilon _{Sr} & \epsilon _{S\mu } & \epsilon _{S\lambda }\\[5pt] 0 & \epsilon _{rr} & \epsilon _{r\mu } & \epsilon _{r\lambda }\\[5pt] 0 & 0 & \epsilon _{\mu \mu } & \epsilon _{\mu \lambda }\\[5pt] 0 & 0 & 0 & 1 \end{array}\right )\quad,\quad \left (\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} 1 & \rho _{Sr} & \rho _{S\mu } & \rho _{S\lambda }\\[5pt] \rho _{Sr} & 1 & \rho _{r\mu } & \rho _{r\lambda }\\[5pt] \rho _{S\mu } & \rho _{r\mu } & 1 & \rho _{\mu \lambda }\\[5pt] \rho _{S\lambda } & \rho _{r\lambda } & \rho _{\mu \lambda } & 1 \end{array}\right )=\Sigma \Sigma ^{\top }\,, \end{eqnarray*}

where $\rho _{Sr,}\,\rho _{S\mu },\,\rho _{S\lambda },\rho _{r\mu },\rho _{r\lambda }$ , and $\rho _{\mu \lambda }\in ({-}1,1)$ . This model incorporates the correlation between financial and mortality shocks, which can be relevant in the context of events such as a pandemic like COVID-19. We assume that costs of risk, noted $\boldsymbol{\theta }=\left (\theta _{j}\right )_{j=1,\ldots,4}^{\top }$ , are constant, i.e., the Brownian motions with drift under $\mathbb{P}$ ,

\begin{eqnarray*} dW_{t}^{(j)} & = & d\tilde{W}_{t}^{(j)}+\theta _{j}dt\,, \end{eqnarray*}

are Brownian motions under $\mathbb{Q}$ . To keep identical dynamics under $\mathbb{P}$ and $\mathbb{Q}$ , the risk premium parameters, stored in a vector $\boldsymbol{\theta }=\left (\theta _{1},\theta _{2},\theta _{3},\theta _{4}\right )^{\top }$ , are such that

\begin{eqnarray*} \nu _{S} & = & \epsilon _{SS}\theta _{1}+\epsilon _{Sr}\theta _{2}+\epsilon _{S\mu }\theta _{3}+\epsilon _{S\lambda }\theta _{4}\,,\\ \nu _{r} & = & -\left (\epsilon _{rr}\theta _{2}+\epsilon _{r\mu }\theta _{3}+\epsilon _{r\lambda }\theta _{4}\right )\,,\\ \nu _{\mu } & = & -\left (\epsilon _{\mu \mu }\theta _{3}+\epsilon _{\mu \lambda }\theta _{4}\right )\,,\\ \nu _{\lambda } & = & -\theta _{4}\,. \end{eqnarray*}

Under the risk neutral measure $\mathbb{Q}$ , financial and biometric processes are ruled by the following SDEs:

(2) \begin{eqnarray} d\left (\begin{array}{c} S_{t}\\ r_{t}\\ \mu _{x_{\mu }+t}\\ \lambda _{x_{\lambda }+t} \end{array}\right ) & = & \left (\begin{array}{c} r_{t}\,S_{t}\\ \kappa _{r}\left (\gamma _{r}(t)-r_{t}\right )\\ \kappa _{\mu }\left (\gamma _{\mu }(t)-\mu _{x_{\mu }+t}\right )\\ \kappa _{\lambda }\left (\gamma _{\lambda }(t)-\lambda _{x_{\lambda }+t}\right ) \end{array}\right )dt+\left (\begin{array}{c} S_{t}\sigma _{S}\Sigma _{S}^{\top }\\[6pt] \sigma _{r}\Sigma _{r}^{\top }\\[6pt] \sigma _{\mu }(t)\Sigma _{\mu }^{\top }\\[6pt] \sigma _{\lambda }(t)\Sigma _{\lambda }^{\top } \end{array}\right )d\boldsymbol{W}_{t}\,. \end{eqnarray}

We will use later these equations to build the valuation equation fulfilled by the GMAB. As mentioned earlier, we have chosen a pure Brownian model specifically for the purpose of benchmarking exact GMAB prices against PINN predictions. The following two sections present intermediate analytical results that are essential for deriving a closed-form expression for the GMAB.

3. Bonds and survival probabilities

We denote by $P(t,T)$ the present value of a zero-coupon bond that pays one monetary unit at expiry date, $T$ . The survival probabilities of an $x_{z}+t$ year old individual up to time $t$ are noted $\,_{T}p_{x_{z}+t}^{z}$ , for $z\in \{\mu,\lambda \}$ . Pure endowments provide a lump sum payment at expiry in case of survival and nothing if the policyholder dies before the term. We denote them by $_{T}E_{t}^{z}$ for $z\in \{\mu,\lambda \}$ . Zero-coupon bonds, survival probabilities, and pure endowents are valued by the following expectations under the risk neutral measure:

\begin{equation*} \begin {cases} P(t,T) & =\mathbb {E}^{\mathbb {Q}}\left (e^{-\int _{t}^{T}r_{s}ds}|\mathcal {F}_{t}\right )\,,\\ \,_{T}p_{x_{z}+t}^{z} & =\mathbb {E}^{\mathbb {Q}}\left (e^{-\int _{t}^{T}z_{x_{z}+s}ds}\,|\,\mathcal {F}_{t}\right )\;z\in \{\mu,\lambda \},\\ _{T}E_{t}^{z} & =N_{t}^{\mu }\mathbb {E}^{\mathbb {Q}}\left (e^{-\int _{t}^{T}(r_{s}+z_{x_{z}+s})ds}\,|\,\mathcal {F}_{t}\right )\;z\in \{\mu,\lambda \}. \end {cases} \end{equation*}

The model being affine, we derive the closed-form expressions of these products. In the remainder of this article, we adopt the following notation

\begin{align*} B_{y}(t,T) & =\frac{1-e^{-y(T-t)}}{y}\,, \end{align*}

where $y\in \mathbb{R}^{+}$ is a positive parameter. Various combined integrals of $B_{y}(t,T)$ , $\sigma _{\mu }(t)$ , and $\sigma _{\lambda }(t)$ are provided in Appendix B and serve to analytically evaluate the GMAB. From Proposition 4 in Appendix A, $\gamma _{r}(T)$ matches the initial yield curve of interest rate if

\begin{eqnarray*} \gamma _{r}(T) & = & -\frac{1}{\kappa _{r}}\partial _{T}^{2}\ln P(0,T)-\partial _{T}\ln P(0,T)+\frac{\sigma _{r}^{2}}{2\kappa _{r}^{2}}\left (1-e^{-2\kappa _{r}T}\right )\,, \end{eqnarray*}

where $-\partial _{T}\ln P(0,T)$ is the instantaneous forward rate. In a similar manner, for a given initial mortality curve $\,_{T}p_{x_{z}}^{z}$ , the function $\gamma _{z}(u)$ is such that

\begin{eqnarray*} \gamma _{z}(T) & = & -\frac{1}{\kappa _{z}}\partial _{T}^{2}\ln \,_{T}p_{x_{z}}^{z}-\partial _{T}\ln \,_{T}p_{x_{z}}^{z}+\frac{\alpha _{z}^{2}e^{2\beta _{z}x_{z}}}{2\kappa _{z}(\kappa _{z}+\beta _{z})}\left (e^{2\beta _{z}T}-e^{-2\kappa _{z}T}\right )\,. \end{eqnarray*}

Bond prices, survival probabilities, and endowments admit closed-form expressions presented in the next proposition.

Proposition 1. The price at time $t\leq T$ of a discount bond of maturity $T$ is linked to the initial interest rate curve at time $t=0$ by the relation

(3) \begin{eqnarray} P(t,T) & = & \exp \left (-r_{t}B_{\kappa _{r}}(t,T)-\left (\partial _{t}\ln P(0,t)\right )B_{\kappa _{r}}(t,T)+\ln \frac{P(0,T)}{P(0,t)}\right )\\ & & \times \exp \left (-\frac{\sigma _{r}^{2}}{4\kappa _{r}}\left (\left (1-e^{-2\kappa _{r}t}\right )B_{\kappa _{r}}(t,T)^{2}\right )\right )\,.\nonumber \end{eqnarray}

In a similar manner, we can show that, when alive at age $x_{z}+t$ , the survival probability up to time $T$ depends on the initial survival term structure as follows for $z\in \{\mu,\lambda \}$

(4) \begin{eqnarray} &&\,_{T}p_{x+t}^{z}=\exp \left (-z_{x_{z}+t}B_{\kappa _{z}}(t,T)-\left (\partial _{t}\ln \,_{t}p_{x_{z}}^{z}\right )B_{\kappa _{z}}(t,T)+\ln \frac{\,_{T}p_{x_{z}}^{z}}{\,_{t}p_{x_{z}}^{z}}\right )\times \\ && \quad \exp \left (\frac{\alpha _{z}^{2}e^{2\beta _{z}(x_{z}+T)}}{2\kappa _{z}^{2}}\left (B_{2\beta _{z}}(t,T)-2B_{2\beta _{z}+\kappa _{z}}(t,T)+B_{2\beta _{z}+2\kappa _{z}}(t,T)\right )\right )\times \nonumber \\ && \quad \exp \left (\frac{\alpha _{z}^{2}e^{2\beta _{z}x_{z}}}{2\kappa _{z}(\kappa _{z}+\beta _{z})}\left (e^{2\beta _{z}T}B_{2\beta _{z}+\kappa _{z}}(t,T)-e^{2\beta _{z}T}B_{2\beta _{z}}(t,T)\right )\right )\times \nonumber \\ && \quad \exp \left (\frac{\alpha _{z}^{2}e^{2\beta _{z}x_{z}}}{2\kappa _{z}(\kappa _{z}+\beta _{z})}\left (e^{-2\kappa _{z}t}B_{2\kappa _{z}}(t,T)-e^{-\kappa _{z}(T+t)}B_{\kappa _{z}}(t,T)\right )\right )\,,\nonumber \end{eqnarray}

whereas the pure endowment contracts are valued by

(5) \begin{eqnarray} && _{T}E_{t}^{\mu }=N_{t}^{\mu }\:_{T}p_{x_{\mu }+t}^{\mu }\:P(t,T)\,\exp \left (\frac{\sigma _{r}\Sigma _{r}^{\top }\Sigma _{\mu }\alpha _{\mu }e^{\beta _{\mu }\left (x_{\mu }+T\right )}}{\kappa _{\mu }\kappa _{r}}\right )\\ & & \quad \times \exp \left (B_{\beta _{\mu }}(t,T)-B_{\kappa _{\mu }+\beta _{\mu }}(t,T)-B_{\kappa _{r}+\beta _{\mu }}(t,T)+B_{\kappa _{\mu }+\kappa _{r}+\beta _{\mu }}(t,T)\right )\,.\nonumber \end{eqnarray}

and

(6) \begin{eqnarray} & & _{T}E_{t}^{\lambda }=N_{t}^{\lambda }\:_{T}p_{x_{\lambda }+t}^{\lambda }\:P(t,T)\,\exp \left (\frac{\sigma _{r}\Sigma _{r}^{\top }\Sigma _{\lambda }\alpha _{\lambda }e^{\beta _{\lambda }\left (x_{\lambda }+T\right )}}{\kappa _{\lambda }\kappa _{r}}\right )\\ & & \quad \times \exp \left (B_{\beta _{\lambda }}(t,T)-B_{\kappa _{\lambda }+\beta _{\lambda }}(t,T)-B_{\kappa _{r}+\beta _{\lambda }}(t,T)+B_{\kappa _{\lambda }+\kappa _{r}+\beta _{\lambda }}(t,T)\right )\,.\nonumber \end{eqnarray}

The sketch of the proof is provided in Appendix A. Using standard developments, the bond price dynamic is equal to

(7) \begin{equation} \frac{dP(t,T)}{P(t,T)}=\left (r_{t}+B_{\kappa _{r}}(t,T)\sigma _{r}\nu _{r}\right )\,dt-B_{\kappa _{r}}(t,T)\sigma _{r}\Sigma _{r}^{\top }d\tilde{\boldsymbol{W}}_{t}\,. \end{equation}

We assume that policyholder’s savings is invested in cash, stocks, zero-coupon bonds, and in mortality bonds written on a reference population of age $x_{\lambda }$ at time $t=0$ . We denote by $D(t,T)$ , the value of this mortality linked bond expiring at time $T$ . Its return depends on the benchmark mortality rate, $\lambda _{x_{\lambda }+t}$ , and is valued at by the following expectation:

\begin{eqnarray*} D(t,T) & = & \mathbb{E^{Q}}\left (e^{-\int _{t}^{T}r_{s}ds-\int _{0}^{T}\lambda _{x_{\lambda }+s}ds}|\mathcal{F}_{t}\right )\\ & = & e^{-\int _{0}^{t}\lambda _{x_{\lambda }+s}ds}\mathbb{E^{Q}}\left (e^{-\int _{t}^{T}\left (r_{s}+\lambda _{x_{\lambda }+s}\right )ds}|\mathcal{F}_{t}\right )\,. \end{eqnarray*}

This product is designed in such a way that, in the event of over- or undermortality, the bondholder realizes a capital loss or gain at maturity. The dynamics of this bond are described in the next proposition, which is proven in Appendix A.

Proposition 2. Under the risk neutral measure $\mathbb{Q}$ , the mortality bond earns on average the risk-free rate and its price obeys to the SDE:

(8) \begin{eqnarray} \frac{dD(t,T)}{D(t,T)} & = & r_{t}dt-\left (B_{\kappa _{r}}(t,T)\sigma _{r}\Sigma _{r}^{\top }+B_{\kappa _{\lambda }}(t,T)\sigma _{\lambda }(t)\Sigma _{\lambda }^{\top }\right )d\boldsymbol{W}_{t} \end{eqnarray}

Under $\mathbb{P}$ , the risk premium of the mortality bond is proportional to interest and mortality volatilities:

(9) \begin{eqnarray} \frac{dD(t,T)}{D(t,T)} & = & \left (r_{t}+B_{\kappa _{r}}(t,T)\sigma _{r}\nu _{r}+B_{\kappa _{\lambda }}(t,T)\sigma _{\lambda }(t)\nu _{\lambda }\right )dt\\ & & -\left (B_{\kappa _{r}}(t,T)\sigma _{r}\Sigma _{r}^{\top }+B_{\kappa _{\lambda }}(t,T)\sigma _{\lambda }(t)\Sigma _{\lambda }^{\top }\right )d\tilde{\boldsymbol{W}}_{t}\nonumber \end{eqnarray}

4. Contract and analytical valuation

GMABs are retirement savings products which promises in case of survival, the maximum between investments, and a guaranteed capital. The policyholder savings are invested in cash, stocks, risk-free and mortality bonds of maturity $T$ . We denote by $\left (A_{t}\right )_{t\geq 0}$ the total value of these investments. The percentages of stocks, risk-free and mortality bonds are assumed constant and denoted by $\boldsymbol{\pi }=\left (\pi _{S},\pi _{P},\pi _{D}\right )_{t\geq 0}$ . We assume that the investment policy is self-financed, the dynamic of assets under management is equal to

\begin{eqnarray*} \frac{dA_{t}}{A_{t}} & = & \pi _{S}\frac{dS_{t}}{S_{t}}+\pi _{P}\frac{dP_{t}}{P_{t}}+\pi _{D}\frac{dD_{t}}{D_{t}}+\left (1-\pi _{S}-\pi _{P}-\pi _{D}\right )r_{t}\,dt\,, \end{eqnarray*}

where $dS_{t}$ , $dP(t,T)$ , and $dD(t,T)$ are, respectively, provided in Equations (1), (7), and (9). We can therefore reformulate the differential of $A_{t}$ under $\mathbb{P}$ as follows

\begin{eqnarray*} \frac{dA_{t}}{A_{t}} & = & \left (r_{t}+\boldsymbol{\pi }^{\top }\boldsymbol{\nu }_{A}(t)\right )dt+\boldsymbol{\pi }^{\top }\Sigma _{A}(t)d\tilde{\boldsymbol{W}}_{t} \end{eqnarray*}

where $\boldsymbol{\nu }_{A}(t)$ is a time-varying vector of risk premiums,

\begin{eqnarray*} \boldsymbol{\nu }_{A}(t) & = & \left (\begin{array}{c} \sigma _{S}\nu _{S}\\ B_{\kappa _{r}}(t,T)\sigma _{r}\nu _{r}\\ B_{\kappa _{r}}(t,T)\sigma _{r}\nu _{r}+B_{\kappa _{\lambda }}(t,T)\sigma _{\lambda }(t)\nu _{\lambda } \end{array}\right )\,, \end{eqnarray*}

and $\Sigma _{A}(t)$ is a volatility $3\times 4$ -matrix :

\begin{equation*} \Sigma _{A}(t)=\left (\begin {array}{c} \sigma _{S}\Sigma _{S}^{\top }\\ -B_{\kappa _{r}}(t,T)\sigma _{r}\Sigma _{r}^{\top }\\ -B_{\kappa _{r}}(t,T)\sigma _{r}\Sigma _{r}^{\top }-B_{\kappa _{\lambda }}(t,T)\sigma _{\lambda }(t)\Sigma _{\lambda }^{\top } \end {array}\right )\,. \end{equation*}

Under the risk neutral measure, investments earn on average the risk-free rate and the differential of $A_{t}$ is ruled by

\begin{eqnarray*} \frac{dA_{t}}{A_{t}} & = & r_{t}dt+\boldsymbol{\pi }^{\top }\Sigma _{A}(t)d\boldsymbol{W}_{t}. \end{eqnarray*}

The contract is subscribed by an individual of age $x_{\mu }$ and guarantees at the expiry date, denoted by $T^{*}$ , a payout equal to the maximum between a guaranteed capital $G_{m}$ and the portfolio $A_{T^{*}}$ , in the event of survival. However, the benefit is eventually bounded by $G_{M}$ . In practice, $G_{m}$ is often set to the initial investment, $A_{0}$ , or to a percentage of $A_{0}$ . Whereas, the participation is most of the time unbounded ( $G_{M}=+\infty$ ). We denote by $\tau _{\mu }\in \mathbb{R}^{+}$ , the random time of insured’s death and $N_{t}^{\mu }=\textbf{1}_{\{\tau _{\mu }\geq t\}}$ . The payoff in case of survival is

(10) \begin{eqnarray} H(A_{T^{*}},G_{m}) & = & \left (G_{T^{*}}+\left (A_{T^{*}}-G_{m}\right )_{+}-\left (A_{T^{*}}-G_{M}\right )_{+}\right ) \end{eqnarray}

The fair value of such a policy, denoted by $L_{t}$ , is equal to the expected discounted cash flows under the chosen risk neutral measure, noted $\mathbb{Q}$ . This is a function of all state variables:

(11) \begin{eqnarray} & & L_{t}=N_{t}^{\mu }\mathbb{E}^{\mathbb{Q}}\left (e^{-\int _{t}^{T^{*}}(r_{s}+\mu _{x_{\mu }+s})ds}G_{m}\,|\,\mathcal{F}_{t}\right )\\ & & \quad \,+N_{t}^{\mu }\mathbb{E}^{\mathbb{Q}}\left (e^{-\int _{t}^{T^{*}}(r_{s}+\mu _{x_{\mu }+s})ds}\left (\left (A_{T^{*}}-G_{m}\right )_{+}\right )\,|\,\mathcal{F}_{t}\right )\nonumber \\ & & \quad \,-N_{t}^{\mu }\mathbb{E}^{\mathbb{Q}}\left (e^{-\int _{t}^{T^{*}}(r_{s}+\mu _{x_{\mu }+s})ds}\left (\left (A_{T^{*}}-G_{M}\right )_{+}\right )\,|\,\mathcal{F}_{t}\right ).\nonumber \end{eqnarray}

In order to obtain a closed form expression of the saving contract (10), we perform a change of measure using as Radon–Nykodym derivative:

(12) \begin{eqnarray} \left .\frac{d\mathbb{F}}{d\mathbb{Q}}\right |_{T^{*}} & =\mathbb{E}^{\mathbb{Q}}\left (\frac{d\mathbb{F}}{d\mathbb{Q}}|\mathcal{F}_{T^{*}}\right )= & \frac{e^{-\int _{0}^{T^{*}}(r_{s}+\mu _{x_{\mu }+s})ds}}{\mathbb{E}^{\mathbb{Q}}\left (e^{-\int _{0}^{T^{*}}(r_{s}+\mu _{x_{\mu }+s})ds}|\mathcal{F}_{0}\right )}\,. \end{eqnarray}

From Equations (A6), this change of measure is rewritten as follows:

\begin{eqnarray*} & & \left .\frac{d\mathbb{F}}{d\mathbb{Q}}\right |_{T^{*}}=\exp \left (-\frac{\sigma _{r}^{2}\epsilon _{rr}^{2}}{2}\int _{0}^{T^{*}}B_{\kappa _{r}}(u,T^{*})^{2}du-\sigma _{r}\epsilon _{rr}\int _{0}^{T^{*}}B_{\kappa _{r}}(u,T^{*})\,dW_{u}^{(2)}\right )\\ & & \quad \times{{\exp \left (-\int _{0}^{T^{*}}\left (\sigma _{r}\epsilon _{r\mu }B_{\kappa _{r}}(u,T^{*})+\sigma _{\mu }(u)\epsilon _{\mu \mu }B_{\kappa _{\mu }}(u,T^{*})\right )dW_{u}^{(3)}\right )}}\\ & & \quad \times \exp \left (-\frac{1}{2}\int _{0}^{T^{*}}\left (\sigma _{r}\epsilon _{r\mu }B_{\kappa _{r}}(u,T^{*})+\sigma _{\mu }(u)\epsilon _{\mu \mu }B_{\kappa _{\mu }}(u,T^{*})\right )^{2}du\right )\\ & & \quad \times{{\exp \left (-\int _{0}^{T^{*}}\left (\sigma _{r}\epsilon _{r\lambda }B_{\kappa _{r}}(u,T^{*})+\sigma _{\mu }(u)\epsilon _{\mu \lambda }B_{\kappa _{\mu }}(u,T^{*})+\sigma _{\lambda }(u)B_{\kappa _{\lambda }}(u,T^{*})\right )dW_{u}^{(4)}\right )}}\\ & & \quad \times{{\exp \left (-\frac{1}{2}\int _{0}^{T^{*}}\left (\sigma _{r}\epsilon _{r\lambda }B_{\kappa _{r}}(u,T^{*})+\sigma _{\mu }(u)\epsilon _{\mu \lambda }B_{\kappa _{\mu }}(u,T^{*})+\sigma _{\lambda }(u)B_{\kappa _{\lambda }}(u,T^{*})\right )^{2}du\right )}} \end{eqnarray*}

We recognize a Doleans–Dade exponential and then under the measure $\mathbb{F}.$ We denote by $\boldsymbol{W}_{t}^{\mathbb{F}}$ the vector of $\mathbb{F}$ -Brownian motions $\left (W_{t}^{(1)\mathbb{F}},W_{t}^{(2)\mathbb{F}},W_{t}^{(3)\mathbb{F}},W_{t}^{(4)\mathbb{F}}\right )$ , defined by

(13) \begin{eqnarray} d\boldsymbol{W}_{t}^{\mathbb{F}} & = & d\boldsymbol{W}_{t}+\left (\sigma _{r}B_{\kappa _{r}}(t,T^{*})\Sigma _{r}+\sigma _{\mu }(t)B_{\kappa _{\mu }}(t,T^{*})\Sigma _{\mu }+\sigma _{\lambda }(t)B_{\kappa _{\lambda }}(t,T^{*})\Sigma _{\lambda }\right )dt\,. \end{eqnarray}

The dynamic of the total asset under the $\mathbb{F}$ - measure is then equal to

(14) \begin{align} & \frac{dA_{t}}{A_{t}}=r_{t}dt+\boldsymbol{\pi }^{\top }\Sigma _{A}(t)d\boldsymbol{W}_{t}^{\mathbb{F}}-\boldsymbol{\pi }^{\top }\Sigma _{A}(t)\times \\ & \left (\sigma _{r}B_{\kappa _{r}}(t,T^{*})\Sigma _{r}+\sigma _{\mu }(t)B_{\kappa _{\mu }}(t,T^{*})\Sigma _{\mu }+\sigma _{\lambda }(t)B_{\kappa _{\lambda }}(t,T^{*})\Sigma _{\lambda }\right )dt\nonumber \end{align}

As $A_{T^{*}}=\frac{A_{T^{*}}}{P(T^{*},T^{*})}$ , we focus on the dynamics of $d\frac{A_{t}}{P(t,T^{*})}$ in order to prove that this ratio is a log-normal process.

Proposition 3. The quantity $\frac{A_{t}}{P(t,T^{*})}$ has a log-normal distribution, $\ln \frac{A_{T^{*}}/P(T^{*},T^{*})}{A_{t}/P(t,T^{*})}$ ∼ $N\left (\mu _{\mathbb{F}}(t),v_{\mathbb{F}}(t)\right )$ , with a time dependent variance equal to

(15) \begin{eqnarray} & & v_{\mathbb{F}}(t)=\int _{t}^{T^{*}}\left (B_{\kappa _{r}}(u,T^{*})\sigma _{r}\Sigma _{r}^{\top }+\boldsymbol{\pi }^{\top }\Sigma _{A}(u)\right )\left (B_{\kappa _{r}}(u,T^{*})\sigma _{r}\Sigma _{r}+\Sigma _{A}(u)^{\top }\boldsymbol{\pi }\right )du\,, \end{eqnarray}

and a time dependent mean given by

(16) \begin{eqnarray} & & \mu _{\mathbb{F}}(t)=-\frac{1}{2}v_{\mathbb{F}}(t)-\sigma _{r}\Sigma _{r}^{\top }\Sigma _{\mu }\int _{t}^{T^{*}}\sigma _{\mu }(u)B_{\kappa _{r}}(u,T^{*})B_{\kappa _{\mu }}(u,T^{*})\,du\\ & & \qquad -\sigma _{r}\Sigma _{r}^{\top }\Sigma _{\lambda }\int _{t}^{T^{*}}\sigma _{\lambda }(u)B_{\kappa _{r}}(u,T^{*})B_{\kappa _{\lambda }}(u,T^{*})du\nonumber \\ & & \qquad -\boldsymbol{\pi }^{\top }\int _{t}^{T^{*}}\Sigma _{A}(u)\sigma _{\mu }(u)B_{\kappa _{\mu }}(u,T^{*})du\,\Sigma _{\mu }\nonumber \\ & & \qquad -\boldsymbol{\pi }^{\top }\int _{t}^{T^{*}}\Sigma _{A}(u)\sigma _{\lambda }(u)B_{\kappa _{\lambda }}(u,T^{*})du\,\Sigma _{\lambda }\,.\nonumber \end{eqnarray}

The proof is provided in Appendix A, while the full analytical expressions of $v_{\mathbb{F}}(t)$ and $\mu _{\mathbb{F}}(t)$ are detailed in Appendix C. As $\frac{A_{t}}{P(t,T^{*})}$ is log-normal under the forward measure, we can evaluate the call options embedded in the GMAB.

Corollary 1. Let us denote by $\Phi (.)$ is the cumulative distribution function (cdf) of a $N(0,1)$ . If we define

(17) \begin{equation} \begin{cases} d_{2}(t)=\frac{\ln \left (\frac{G_{T^{*}}}{A_{t}/P(t,T^{*})}\right )-\mu _{\mathbb{F}}(t)}{\sqrt{v_{\mathbb{F}}(t)}} & \,,\,d_{1}(t)=d_{2}(t)-\sqrt{v_{\mathbb{F}}(t)}\negmedspace,\\ d_{2}^{M}(t)=\frac{\ln \left (\frac{G_{M}}{A_{t}/P(t,T^{*})}\right )-\mu _{\mathbb{F}}(t)}{\sqrt{v_{\mathbb{F}}(t)}} & \,,\,d_{1}^{M}(t)=d_{2}^{M}(t)-\sqrt{v_{\mathbb{F}}(t)}\negmedspace \end{cases} \end{equation}

then expected positive differences between $A_{T^{*}}$ , $G_{T^{*}}$ , and $G_{M}$ under the forward measure are given by

\begin{equation*} \begin {array}{c} \mathbb {E}^{\mathbb {F}}\left (\left (A_{T^{*}}-G_{m}\right )_{+}|\mathcal {F}_{t}\right )=\frac {A_{t}}{P(t,T^{*})}e^{\mu _{\mathbb {F}}(t)+\frac {v_{\mathbb {F}}(t)}{2}}\Phi \left (-d_{1}(t)\right )-G_{m}\Phi \left (-d_{2}(t)\right )\,,\\ \mathbb {E}^{\mathbb {F}}\left (\left (A_{T^{*}}-G_{M}\right )_{+}|\mathcal {F}_{t}\right )=\frac {A_{t}}{P(t,T^{*})}e^{\mu _{\mathbb {F}}(t)+\frac {v_{\mathbb {F}}(t)}{2}}\Phi \left (-d_{1}^{M}(t)\right )-G_{M}\Phi \left (-d_{2}^{M}(t)\right )\,. \end {array} \end{equation*}

This corollary is proven using standard Black & Scholes developments. This last result allows us to infer that the market value of the GMAB (11) is given by

(18) \begin{eqnarray} & & L_{t}=\,_{T^{*}}E_{t}^{\mu }\,G_{m}\\ & & +\,_{T^{*}}E_{t}^{\mu }\left [\frac{A_{t}e^{\mu _{\mathbb{F}}(t)+\frac{v_{\mathbb{F}}(t)}{2}}}{P(t,T^{*})}\Phi \left (-d_{1}(t)\right )-G_{m}\,\Phi \left (-d_{2}(t)\right )\right ]\nonumber \\ & & -\,_{T^{*}}E_{t}^{\mu }\left [\frac{A_{t}e^{\mu _{\mathbb{F}}(t)+\frac{v_{\mathbb{F}}(t)}{2}}}{P(t,T^{*})}\Phi \left (-d_{1}^{M}(t)\right )-G_{M}\Phi \left (-d_{2}^{M}(t)\right )\right ]\,,\nonumber \end{eqnarray}

where $\,_{T^{*}}E_{t}^{\mu }$ is given by Equation (5). This closed form is used to measure the accuracy of the PINN for predicting GMAB prices.

5. Scaled Feynman–Kac equation

A PINN integrates principles from physics, including PDEs governing the behavior of state variables. In a financial context, we adopt the Feynman–Kac (FK) equation as our guiding principle. The FK equation is a PDE satisfied by all assets traded in an arbitrage-free market. In this section, we construct this equation for the model under consideration in this article. In the following section, we utilize a neural network to solve it for various investment strategies, contract maturities, and bond maturities. To lighten future developments, we denote the vector of state variables, augmented with the total asset value, by

\begin{equation*} \boldsymbol{y}_{t}=\left (S_{t},r_{t},\mu _{x_{\mu }+t},\lambda _{x_{\lambda }+t},A_{t}\right )^{\top }\,. \end{equation*}

Under the risk neutral measure $\mathbb{Q}$ , this multivariate process is ruled by the following SDE:

\begin{eqnarray*} \text{d}\boldsymbol{y}_{t} & = & \boldsymbol{\mu} _{\boldsymbol{y}}(t,\boldsymbol{y}_{t})dt+\Sigma _{\boldsymbol{y}}(t,\boldsymbol{y}_{t})d\boldsymbol{W}_{t}\,, \end{eqnarray*}

where $\boldsymbol{\mu} _{\boldsymbol{y}}(.)$ is a vector of dimension 5 and $\Sigma _{\boldsymbol{y}}(.)$ is a $5\times 4$ matrix:

\begin{equation*} \boldsymbol{\mu}_{\tilde{\boldsymbol{y}}}(t,\boldsymbol{y}_{t})=\left (\begin {array}{c} \begin {array}{c} \begin {array}{c} r_{t}\,S_{t}\\ \kappa _{r}\left (\gamma _{r}(t)-r_{t}\right )\\ \kappa _{\mu }\left (\gamma _{\mu }(t)-\mu _{x_{\mu }+t}\right )\\ \kappa _{\lambda }\left (\gamma _{\lambda }(t)-\lambda _{x_{\lambda }+t}\right )\\ r_{t}A_{t} \end {array}\end {array}\end {array}\right )\;,\;\Sigma _{\boldsymbol{y}}(t,\boldsymbol{y}_{t})=\left (\begin {array}{c} S_{t}\sigma _{S}\Sigma _{S}^{\top }\\[5pt] \sigma _{r}\Sigma _{r}^{\top }\\[5pt] \sigma _{\mu }(t)\Sigma _{\mu }^{\top }\\[5pt] \sigma _{\lambda }(t)\Sigma _{\lambda }^{\top }\\[5pt] A_{t}\boldsymbol {\pi }^{\top }\Sigma _{A}(t) \end {array}\right )\,. \end{equation*}

The GMAB is a function of time and state variables parameterized by the investment policy $\boldsymbol{\pi }$ , the duration of the contract $T^{*}$ , and the maturity of bonds, $T$ . As we aim to understand the relation between assets and liabilities, we emphasize the dependence to investment parameters by denoting the contract price as follows:

\begin{equation*} L_{t}=V\left (t,\boldsymbol{y}_{t},N_{t}^{\mu }\,|\,\boldsymbol {\pi },T^{*},T\right )\,. \end{equation*}

Under the assumption of arbitrage-free market, all traded securities, including the insurance contract, earn on average the risk-free rate, i.e. :

(19) \begin{align} \mathbb{E}\left (dL_{t+}\,|\,\mathcal{F}_{t}\right ) & =L_{t}\,r_{t}\,dt\,. \end{align}

Let us, respectively, denote the gradient and the Hessian of $V(.)$ with respect to $\boldsymbol{y}$ by $\nabla _{\boldsymbol{y}}V$ and $\mathcal{H}_{\boldsymbol{y}}(V)$ :

\begin{eqnarray*} \nabla _{\boldsymbol{y}}V=\left (\begin{array}{c} \partial _{S}V\\ \partial _{r}V\\ \partial _{\mu }V\\ \partial _{\lambda }V\\ \partial _{A}V \end{array}\right ) & \:,\: & \mathcal{H}_{\boldsymbol{y}}(V)=\left (\begin{array}{ccccc} \partial _{SS}V & \partial _{Sr}V & \partial _{S\mu }V & \partial _{S\lambda }V & \partial _{SA}V\\ \partial _{Sr}V & \partial _{rr}V & \partial _{r\mu }V & \partial _{r\lambda }V & \partial _{rA}V\\ \partial _{S\mu }V & \partial _{r\mu }V & \partial _{\mu \mu }V & \partial _{\mu \lambda }V & \partial _{\mu A}V\\ \partial _{S\lambda }V & \partial _{r\lambda }V & \partial _{\mu \lambda }V & \partial _{\lambda \lambda }V & \partial _{\lambda A}V\\ \partial _{SA}V & \partial _{rA}V & \partial _{\mu A}V & \partial _{\lambda A}V & \partial _{AA}V \end{array}\right )\,. \end{eqnarray*}

Applying the Itô’s lemma to $V\left (.\right )$ allows us to rewrite (19) as a partial differential valuation equation:

(20) \begin{align} 0= & \partial _{t}V-\left (r_{t}+\mu _{x+t}\right )\,V+\boldsymbol{\mu} _{\boldsymbol{y}}(t,\boldsymbol{y}_{t})^{\top }\nabla _{\boldsymbol{y}}V\\ & +\frac{1}{2}\text{tr}\left (\Sigma _{\boldsymbol{y}}(t,\boldsymbol{y}_{t})\Sigma _{\boldsymbol{y}}(t,\boldsymbol{y}_{t})^{\top }\mathcal{H}_{\boldsymbol{y}}(V)\right )\,.\nonumber \end{align}

This last expression is called the “Feynman–Kac (FK) equation”. The boundary constraints on $V(.)$ are

(21) \begin{equation} \begin{cases} V(T^{*},\boldsymbol{y}_{T^{*}},1\,|\,\boldsymbol{\pi },T^{*},T) & =H\left (A_{T^{*}},G_{m},G_{M}\right )\,,\\ V(t,\boldsymbol{y}_{t},0\,|\,\boldsymbol{\pi },T^{*},T) & =0\,. \end{cases} \end{equation}

Notice that we have 5 state variables and 4 Brownian motions. The model is, therefore, over-specified, and we can remove one state variable from $\boldsymbol{y}$ , for example, the stock price, which is redundant with respect to the total asset value, $A_{t}$ . Nonetheless, we have retained all state variables in the numerical illustration to observe how the PINN handles this naive over-specification.

Due to the numerous risk factors involved, solving numerically the PDE (20) is numerically challenging. Implementing either an explicit or implicit finite difference method necessitates constructing a 6-dimensional grid in the spacetime of time and state variables, with relatively fine meshes to ensure convergence. In practice, this becomes intractable. As an alternative, we approximate the solution to the FK Equation (20) using a neural network that takes as input the time, state variables, asset mix, and durations of GMAB and bonds. The mathematical definition of a neural network will be revisited in the following section, but it is worth noting that a neural network employs bounded activation functions, such as sigmoid or hyperbolic tangents. To mitigate convergence issues associated with the vanishing gradient problem, we must standardize and scale the network’s inputs. This preprocessing step slightly modifies the FK equation. Let us define the vectors $\boldsymbol{a},\boldsymbol{b}\in \mathbb{R}^{5}$ as follows

(22) \begin{eqnarray} \boldsymbol{a} & = & \left (a_{S},a_{r},a_{\mu },a_{\lambda },a_{A}\right )^{\top }\,,\\ \boldsymbol{b} & = & \left (b_{S},b_{r},b_{\mu },b_{\lambda },b_{A}\right )^{\top }\,.\nonumber \end{eqnarray}

$\boldsymbol{a}$ and $\boldsymbol{b}$ are chosen to normalize a random sample of state variables. This point is discussed in the next section. The vector of centered and scaled state variables is denoted by

(23) \begin{eqnarray} \tilde{\boldsymbol{y}_{t}} & = & \boldsymbol{a}+\boldsymbol{b}\odot \boldsymbol{y}_{t}\,, \end{eqnarray}

where $\odot$ is the elementwise product. In a similar manner, we center and scale the time with coefficients $a_{h}$ and $b_{h}$ that will depends on the horizon of valuation. The centered scaled time and durations are defined as follows:

(24) \begin{equation} \begin{cases} \tilde{t} & =a_{h}+b_{h}t\,,\\ \tilde{T} & =a_{h}+b_{h}T\,,\\ \tilde{T}^{*} & =a_{h}+b_{h}T^{*}. \end{cases} \end{equation}

We do not need to scale the vector $\boldsymbol{\pi }\in [0,1]^{3}$ . The normalized state variables are ruled by the SDE:

\begin{eqnarray*} d{{\tilde{\boldsymbol{y}}_{t}}} & = & \boldsymbol{\mu} _{\tilde{\boldsymbol{y}}}(t,{\tilde{\boldsymbol{y}}}_{t})dt+\Sigma _{{\tilde{\boldsymbol{y}}}}(t,{\tilde{\boldsymbol{y}}}_{t})d\boldsymbol{W}_{t}\,, \end{eqnarray*}

where $\mu _{\tilde{y}}(.)$ is a $5\times$ vector and $\Sigma _{{\tilde{\boldsymbol{y}}}}(.)$ is a 5 $\times$ 4 matrix:

\begin{equation*} \boldsymbol{\mu} _{\tilde{\boldsymbol{y}}}(t,\tilde{\boldsymbol{y}}_{t})=\left (\begin{array}{c} \left (\tilde {S}_{t}-a_{S}\right )\frac {\tilde {r}_{t}-a_{r}}{b_{r}}\\ \kappa _{r}\left (b_{r}\gamma _{r}(t)+a_{r}-\tilde {r}_{t}\right )\\ \kappa _{\mu }\left (b_{\mu }\gamma _{\mu }(t)+a_{\mu }-\tilde {\mu }_{x_{\mu }+t}\right )\\ \kappa _{\lambda }\left (b_{\lambda }\gamma _{\lambda }(t)+a_{\lambda }-\tilde {\lambda }_{x_{\lambda }+t}\right )\\ \left (\tilde {A}_{t}-a_{A}\right )\frac {\tilde {r}_{t}-a_{r}}{b_{r}} \end{array}\right )\;,\;\Sigma _{\tilde{\boldsymbol{y}}}(t,\tilde{\boldsymbol{y}}_{t})=\left(\begin{array}{c} \left (\tilde {S}_{t}-a_{S}\right )\sigma _{S}\Sigma _{S}^{\top }\\ b_{r}\sigma _{r}\Sigma _{r}^{\top }\\ b_{\mu }\sigma _{\mu }(t)\Sigma _{\mu }^{\top }\\ b_{\lambda }\sigma _{\lambda }(t)\Sigma _{\lambda }^{\top }\\ \left (\tilde {A}_{t}-a_{A}\right )\boldsymbol {\pi }_{t}^{\top }\Sigma _{A}(t) \end{array}\right)\,. \end{equation*}

The value of the contract may be seen as a function of time and rescaled state variables, $L_{t}=V(\tilde{t},\tilde{\boldsymbol{y}}_{t},N_{t}^{\mu },|\,\boldsymbol{\pi },\tilde{T}^{*},\tilde{T})$ . The valuation Equation (20) is then rewritten in rescaled form:

(25) \begin{align} 0= & b_{h}\,\partial _{\tilde{t}}V-\left (\frac{\tilde{r}_{t}-a_{r}}{b_{r}}+\frac{\tilde{\mu }_{x+t}-a_{\mu }}{b_{\mu }}\right )\,V+\boldsymbol{\mu} _{\tilde{\boldsymbol{y}}}(t,\tilde{\boldsymbol{y}}_{t})^{\top }\nabla _{{\tilde{\boldsymbol{y}}}}V\\ & +\frac{1}{2}\text{tr}\left (\Sigma _{{\tilde{\boldsymbol{y}}}}(t,{\tilde{\boldsymbol{y}}}_{t})\Sigma _{{\tilde{\boldsymbol{y}}}}(t,{\tilde{\boldsymbol{y}}}_{t})^{\top }\mathcal{H}_{{\tilde{\boldsymbol{y}}}}(V)\right )\,,\nonumber \end{align}

where $\nabla _{{\tilde{\boldsymbol{y}}}}V$ and $\mathcal{H}_{{\tilde{\boldsymbol{y}}}}(V)$ are, respectively, the gradient and the Hessian of $V$ with respect to standardized state variables, ${\tilde{\boldsymbol{y}}}$ . The rescaled version of boundary conditions (21) is

(26) \begin{equation} \begin{cases} V\left (\tilde{T}^{*},{\tilde{\boldsymbol{y}}}_{T^{*}},1\,|\,\boldsymbol{\pi },\tilde{T}^{*},\tilde{T}\right ) & =H\left (\left (\tilde{A}_{T^{*}}-a_{A}\right )/b_{A},G_{T^{*}}\right )\,,\\ V\left (\tilde{t},{\tilde{\boldsymbol{y}}}_{t},0\,|\,\boldsymbol{\pi },\tilde{T}^{*},\tilde{T}\right ) & =0\,. \end{cases} \end{equation}

The next section explains how to solve Equation (25) with a neural network. Notice that the stock and total asset values exhibit geometric dynamics. An alternative formulation of the FK equation may be derived by considering log-stock and log-asset prices as state variables. This has the advantage of partially reducing the nonlinearity of the FK PDE. However, in the existing literature, such a transformation is not applied and does not adversely affect the accuracy of PINNs. For instance, Sirignano & Spiliopoulos (Reference Sirignano and Spiliopoulos2018), Glau & Wunderlich (Reference Glau and Wunderlich2022), or Gatta et al. (Reference Gatta, Di Cola, Giampaolo and Piccialli F.Cuomo2023) evaluate European or American basket options by solving the Feynman–Kac equation in a multivariate geometric market without prior log-transform.

In practice, a log-transform is deemed unnecessary because the partial derivatives of the approximated GMAB prices are computed precisely through automatic differentiation, not finite differences. This ensures the robustness of the valuation of the Feynman–Kac equation.

Figure 1 Feed forward network with skip connections toward intermediate layers.

6. Neural networks

We approximate the value function solving Equation (25) using a particular type of neural network. This network takes as input a vector of dimension 11, denoted by ${\boldsymbol{z}}$ , containing the scaled time and state variables, the asset mix vector, and scaled contract and bond maturities:

(27) \begin{eqnarray} {\boldsymbol{z}} & :\!= & \left (\tilde{t},{\tilde{\boldsymbol{y}}}_{t},\boldsymbol{\pi },\tilde{T}^{*},\tilde{T}\right )^{\top }\,. \end{eqnarray}

There is no systematic procedure for determining the optimal structure of a neural network. As proven by Hornik (Reference Hornik1991), we know that single-layer neural networks can approximate regular functions arbitrarily well, but achieving reasonable accuracy may require a high number of neurons. An alternative approach is provided by feed-forward networks, where information flows forward through several neural layers. However, these networks often struggle to replicate non-linear functions with a limited number of layers.

To address this challenge, we adopt a network with an architecture illustrated in Fig. 1 This network is a basic version of recurrent neural networks. During the valuation of the network response, the outputs of neurons in hidden layers, added to the initial data vector, serve as inputs for the immediate next layer. The long short-term memory (LSTM) networks are built on the same principle and are powerful tools for time-series prediction (see, e.g., Denuit et al., Reference Denuit, Hainaut and Trufin2019, Chapter 8). The deep Galerkin method (DGM) network, initially proposed by Sirignano & Spiliopoulos (Reference Sirignano and Spiliopoulos2018) for PINN pricing of basket options, belongs to the family of LSTM models. Recurrent neural networks better replicate non-linear functions than classical feed-forward networks. The same network as the one proposed in this article is successfully applied in Hainaut & Casas (Reference Hainaut and Casas2024) for pricing European options in the Heston model.

The network used in this work presents the same drawbacks as a classical recurrent architecture. A model with too many layers becomes untrainable due to the phenomenon of the “vanishing gradient,” as revealed by Hochreiter et al. (Reference Hochreiter, Bengio, Frasconi and Schmidhuber2001). The lagged responses from hidden neural layers pass through several intermediate activation functions, having an indirect feedback effect on the final response of the network. This effect may become marginal when considering too many intermediate layers, rendering the network untrainable.

Definition Let $l$ , $n_{0}$ , $n_{1}$ , …, $n_{l}\in \mathbb{N}$ be, respectively, the number of layers and neurons in each layer. $n_{0}$ is also the size of input vector. The activation function of layer $k=1,2,\ldots,l$ is noted $\phi _{k}(.)\,{:}\,\mathbb{R}\to \mathbb{R}$ . Let $C_{1}\in \mathbb{R}^{n_{1}}\times \mathbb{R}^{n_{0}}$ , ${\boldsymbol{c}}_{1}\in \mathbb{R}^{n_{1}}$ , $C_{k}\in \mathbb{R}^{n_{k}}\times \mathbb{R}^{n_{0}+n_{k-1}}$ , ${\boldsymbol{c}}_{k}\in \mathbb{R}^{n_{k}}$ for $k=2,\ldots,l-1$ , $C_{l}\in \mathbb{R}^{n_{l}}\times \mathbb{R}^{n_{l-1}}$ , ${\boldsymbol{c}}_{l}\in \mathbb{R}^{n_{l}}$ be neural weights defining the input, intermediate, and output layers. We define the following functions

\begin{equation*} \begin {cases} \psi _{k}({\boldsymbol{x}})=\phi _{k}\left (C_{k}{\boldsymbol{x}}+{\boldsymbol{c}}_{k}\right )\,, & k=1,l\\ \psi _{k}({\boldsymbol{x}},{\boldsymbol{z}})=\phi _{k}\left (C_{k}\left (\begin {array}{c} {\boldsymbol{x}}\\ {\boldsymbol{z}} \end {array}\right )+{\boldsymbol{c}}_{k}\right )\,, & k=2,\ldots,l-1 \end {cases} \end{equation*}

where activation functions $\phi _{k}(.)$ are applied componentwise. The residual neural network is a function $F\,{:}\,\mathbb{R}^{n_{0}}\to \mathbb{R}^{n_{l}}$ defined by

\begin{eqnarray*} F({\boldsymbol{z}}) & = & \psi _{l-1}\circ \psi _{l-1}\left (.,{\boldsymbol{z}}\right )\ldots \circ \psi _{2}\left (.,{\boldsymbol{z}}\right )\circ \psi _{1}({\boldsymbol{z}})\,. \end{eqnarray*}

After having chosen a network architecture, the model is trained by minimizing a loss function that is proportional to errors of approximation. This error is measured by replacing $V(.)$ with $F(.)$ in the scaled FK Equation (25), at random points in the domain.

At time $t\in [0,T^{*}]$ , the domain of the state vector $\boldsymbol{y}_{t}=\left (S_{t},r_{t},\mu _{x_{\mu }+t},\lambda _{x_{\mu }+t},A_{t}\right )^{\top }$ is $\,\mathbb{R}^{+}\times \mathbb{R}\times \mathbb{R}^{+\,3}$ . We approximate this domain by a closed convex subspace:

\begin{align*} \mathcal{D}^{\boldsymbol{y}}\, & \,{:}\,[S_{l},S_{u}]\times [r_{l},r_{u}]\times [\mu _{l},\mu _{u}]\times [\lambda _{l},\lambda _{u}]\times [A_{l},A_{u}]\,. \end{align*}

In order to fit the neural network, we draw a sample of $n_{\mathcal{D}}$ realizations of $\boldsymbol{y}_{t}$ in $\mathcal{D}^{\boldsymbol{y}}$ , at random times. We also sample parameters $\left (\boldsymbol{\pi }_{j},T_{j}^{*},T_{j}\right )_{j=1,\ldots,,n_{\mathcal{D}}}$ under the constraints:

\begin{equation*} \begin {cases} \boldsymbol {\pi }_{j}\in [0,1]\times [0,1]\times [0,1]\,,\\ T_{j}^{*}\in [0,T_{max}^{*}]\,,\,T_{j}\in [0,T_{max}]\\ t_{j}\leq T_{j}^{*}\leq T_{j}\,. \end {cases} \end{equation*}

The set of sampled state variables and parameters is noted $\mathcal{S}_{\mathcal{D}}=\left (t_{j},\boldsymbol{y}_{j},\pi _{j},T_{j}^{*},T_{j}\right )_{j=1,\ldots,n_{\mathcal{D}}}$ . We next center and scale state variables and times. The mean and standard deviation of sampled state variables are computed with the standard formulas:

\begin{equation*} \bar {\boldsymbol{y}}=\frac {1}{n_{\mathcal {D}}}\sum _{j=1}^{n_{\mathcal {D}}}\boldsymbol{y}_{j}\quad,\quad {\boldsymbol{S}}_{y}=\sqrt {\frac {1}{n_{\mathcal {D}}-1}\sum _{j=1}^{n_{\mathcal {D}}}\left (\boldsymbol{y}_{j}-\bar{\boldsymbol{y}}\right )^{2}}\,. \end{equation*}

The scaling vectors (22) are defined by $\boldsymbol{a}=-\frac{\bar{\boldsymbol{y}}}{{\boldsymbol{S}}_{\boldsymbol{y}}}$ and $\boldsymbol{b}=\frac{1}{{\boldsymbol{S}}_{\boldsymbol{y}}}$ . For $1,\ldots,n_{\mathcal{D}}$ , we define $\tilde{\boldsymbol{y}_{j}}=\boldsymbol{a}+\boldsymbol{b}\odot \boldsymbol{y}_{j}\,$ . The times, contract, and bond maturities are also rescaled with relations (24) and coefficients

\begin{eqnarray*} a_{h}= & -\frac{1}{2}\:,\: & b_{h}=\frac{1}{T_{max}^{*}}\,. \end{eqnarray*}

The set of sampled normalized state variables and parameters is denoted by

\begin{equation*} \tilde {\mathcal {S}}_{\mathcal {D}}=\left (\tilde {t}_{j},\tilde {\boldsymbol{y}}_{j},\pi _{j},\tilde {T}_{j}^{*},\tilde {T}_{j}\right )_{j=1,\ldots,n_{\mathcal {D}}}\,. \end{equation*}

During the training phase, the error of approximation is measured for all points of this sample sets with the scaled FK Equation (26). The error at expiry is measured with another sample set, denoted by $\mathcal{\tilde{S}}_{T^{*}}$ , of $n_{T^{*}}$ realizations of state variables and parameters:

\begin{equation*} \mathcal {\tilde {S}}_{T^{*}}=\left (\tilde {\boldsymbol{y}}_{k},\pi _{k},\tilde {T}_{k}^{*},\tilde {T}_{k}\right )_{k=1,\ldots,n_{T^{*}}}\,. \end{equation*}

Let us denote by $\boldsymbol{\Theta },$ the vector containing all neural weights $\left (C_{k},{\boldsymbol{c}}_{k}\right )_{k=1,\ldots,l}$ . At points of $\tilde{\mathcal{S}}_{\mathcal{D}}$ , we define for $j=1,\ldots,n_{\mathcal{D}}$ , the error in Equation (25) when $V$ is replaced by the neural network as follows:

(28) \begin{align} & e_{j}^{\mathcal{D}}(\boldsymbol{\Theta })=b_{h}\,\partial _{\tilde{t}_{j}}F-\left (\frac{\tilde{r}_{j}-a_{r}}{b_{r}}+\frac{\tilde{\mu }_{j}-a_{\mu }}{b_{\mu }}\right )\,F\\ & \quad +\boldsymbol{\mu} _{\tilde{\boldsymbol{y}}}(t_{j},\tilde{\boldsymbol{y}}_{j})^{\top }\nabla _{{\tilde{\boldsymbol{y}}}}F+\frac{1}{2}\text{tr}\left (\Sigma _{{\tilde{\boldsymbol{y}}}}(t_{j},\tilde{\boldsymbol{y}}_{j})\Sigma _{{\tilde{\boldsymbol{y}}}}(t_{j},\tilde{\boldsymbol{y}}_{j})^{\top }\mathcal{H}_{{\tilde{\boldsymbol{y}}}}(F)\right )\,.\nonumber \end{align}

We refer to $e_{j}^{\mathcal{D}}$ as the error on the interior domain since it measures the goodness of fit before expiry. The average quadratic loss on $\tilde{\mathcal{S}}_{\mathcal{D}}$ is the first component of the total loss function used to fit the neural network,

(29) \begin{equation} \mathcal{L}_{D}\left (\boldsymbol{\Theta }\right )=\frac{1}{n_{\mathcal{D}}}\sum _{j=1}^{n_{\mathcal{D}}}e_{j}^{\mathcal{D}}(\boldsymbol{\Theta })^{2}\,. \end{equation}

Since the error $e_{j}^{\mathcal{D}}(\boldsymbol{\Theta })$ depends on first- and second-order partial derivatives, special attention must be given to the accuracy of their calculation. Computing these derivatives using a standard finite difference method may introduce numerical instabilities. Therefore, we opt for an alternative approach known as automatic or algorithmic differentiation. Automatic differentiation leverages the fact that every computer calculation executes a sequence of elementary arithmetic operations and elementary functions, allowing us to compute partial derivatives with accuracy up to the working precision. Automatic differentiation is implemented in TensorFlow through functions such as GradientTape(.) and tape.gradient(.). For a more in-depth introduction and perspectives, we refer the reader to van Merriënboer et al. (Reference van Merriënboer, Breuleux and Bergeron2018).

On the boundary sample set $\mathcal{\tilde{S}}_{T^{*}}$ , we define the error $e_{k}^{T^{*}}(\boldsymbol{\Theta })$ for $k=1,\ldots,n_{T^{*}}$ , as the difference between the output of the neural network and the payoff at expiry:

\begin{eqnarray*} e_{k}^{T^{*}}(\boldsymbol{\Theta }) & = & F(\tilde{T}_{k}^{*},{\tilde{\boldsymbol{y}}}_{k},\boldsymbol{\pi }_{k},\tilde{T}_{k}^{*},\tilde{T}_{k})-H\left (\left (\tilde{A}_{k}^{*}-a_{S}\right )/b_{S},G_{k}\right )\,. \end{eqnarray*}

The average quadratic loss at expiry is the second component of the total loss function.

(30) \begin{equation} \mathcal{L}_{T^{*}}\left (\boldsymbol{\Theta }\right )=\frac{1}{n_{T^{*}}}\sum _{k=1}^{n_{T^{*}}}e_{k}^{T^{*}}(\boldsymbol{\Theta })^{2}\,. \end{equation}

The optimal network weights are found by minimizing the losses in $\tilde{\mathcal{S}}_{\mathcal{D}}$ and $\mathcal{\tilde{S}}_{T^{*}}$ ,

(31) \begin{align} \boldsymbol{\Theta }_{opt} & =\arg \min _{\boldsymbol{\Theta }}\left [\mathcal{L}_{D}\left (\boldsymbol{\Theta }\right )+\mathcal{L}_{T^{*}}\left (\boldsymbol{\Theta }\right )\right ]\:. \end{align}

In practice, the minimization is performed with a gradient descent algorithm. For an introduction, we refer, e.g., to Denuit et al. (Reference Denuit, Hainaut and Trufin2019), Section 1.6.

7. Numerical illustration

We fit a Nelson–Siegel model to the Belgian state yield curveFootnote 1 on the 1/9/23. Initial survival probabilities, $\,_{t}p_{x}^{\mu }$ and $\sigma _{\mu }(t)$ , are fitted to male Belgian mortality rates.Footnote 2 Details are provided in Appendix D. Parameters of survival probabilities $\,_{t}p_{x}^{\lambda }$ , of the reference population for the mortality bond, are assumed proportional to these of $\,_{t}p_{x}^{\mu }$ . Reference ages, $x_{\mu }$ and $x_{\lambda }$ , are set to $50$ years. Other market parameters are estimated from daily time series of the Belgian stock index BEL 20 and of the 1 year Belgian state yield from the 1/9/10 to the 1/9/23. The correlations $\rho _{S\,\mu }$ and $\rho _{r\,\mu }$ are set to –5% and 0%. Parameter estimates are reported in Table 1.

Table 1. Financial and biometric parameters

We first build a sample set, $\mathcal{S_{D}}$ , of contract features, asset-mixes, and durations. We randomly draw $n_{\mathcal{D}}=100,000$ combinations of risk factors and investment-duration parameters. Most of these quantities are drawn from uniform distributions whose characteristics are reported in Table 2. The contract and bond maturities are drawn from exponential distributions to ensure a smoother allocation of durations than what would be obtained with uniform distributions. As $\pi _{S}$ , $\pi _{P}$ , and $\pi _{D}$ are randomly distributed in $[0,1]$ , we allow for short positions in cash. Using the same distributions as for $\mathcal{S_{D}}$ , we generate a boundary sample set $\mathcal{S}_{T^{*}}$ of size, $n_{T^{*}}=50,000$ . The last step consists in scaling and centering the sample sets, as described in the previous section. Notice that setting $G_{m}=100$ is not constraining. As the payoff is piecewise linear, we can use the rule of thumb to evaluate any contract with another minimum capital. If $V(A_{t},G_{m})$ is the GMAB price at time $t$ , for an asset value $A_{t}$ and a guarantee $G_{m}$ , the value of a contract with a guaranteed capital $G_{m}^{'}\neq G_{m}$ is equal to:

Table 2. Model parameters and features of contracts

\begin{equation*} V(A_{t},G_{m}^{'})=\frac {G_{m}^{'}}{G_{m}}V\left (\frac {G_{m}}{G_{m}^{'}}A_{t},G_{m}\right )\,. \end{equation*}

The neural network is implemented in TensorFlow using Python. As mentioned earlier, the necessary partial derivatives for evaluating the inner loss are computed through automatic differentiation. This procedure is implemented in TensorFlow functions like GradientTape(.) and tape.gradient(.), allowing for accurate derivative calculations without relying on finite differences. To train the network, we employ the ADAM optimizer and utilize random batches containing 34,000 and 17,000 samples from $\mathcal{\tilde{S}_{D}}$ and $\mathcal{\tilde{S}}_{T^{*}}$ , respectively. Over the course of one epoch, neural weights are updated three times. We run a total of 26,000 epochs and adjust the ADAM learning rate based on the pattern provided in Table 3, to expedite calibration. The training is conducted on a laptop equipped with an AMD Ryzen 7 5825U processor and 16 Gb of RAM (without a CUDA compatible GPU). The training time depends on the network configuration. As indicated in Table 4, running 200 epochs takes approximately 12–18 minutes for networks with 3 and 5 hidden layers with 32 neurons, respectively. The entire calibration process spans from 26 to 40 hours. Choosing 32 neurons per hidden layer strikes a good balance between accuracy and training time. Tests conducted with 16 neurons over 8000 epochs revealed that this configuration performed less effectively compared to the 32-neuron networks.

Table 3. Learning rate in function of Epochs

Table 4. Information about training time and loss values after training for three neural configurations

Table 4 reports the losses on $\mathcal{\tilde{S}_{D}}$ , $\mathcal{\tilde{S}}_{T^{*}}$ and their total, computed with Equations (29) and (30). We observe a significant reduction of losses between 3 and 4 hidden layers models with 32 neurons. The reduction is relatively smaller when we switch from a 4 to a 5 layers model.

Table 5 presents statistics about relative errorsFootnote 3 between PINN and exact prices obtained with formula (18). These statistics are computed for the 100,000 contracts in the training data set $\tilde{\mathcal{S}}_{\mathcal{D}}$ and for a validation set of same size and built in the same manner. The similarity of figures on both sets does not reveal that we overfit the training dataset. The average relative errors are small and around 0.15% for 4 and 5 hidden layers models. The 95% confidence level is also below 1% for these networks which is quite acceptable.

Table 5. Relative errors computed on the 100 000 contracts of the training set $\tilde{\mathcal{S}}_{\mathcal{D}}$ and on a validation set of same size

In the remainder of this section, we focus on the network with 4 hidden layers as it yields the lowest relative errors on training and validation sets. Fig. 2 compares exact and approximated PINN prices for an 8-year GMAB. The upper plots show the values at issuance and at expiry for various values of the total asset. Market conditions are as of September 1, 2023. Across a wide range of asset values, the exact and PINN prices closely align. The largest deviations are observed for the lowest and highest values of $A_{t}$ . The lower plots display the value at issuance of the same contract with $A_{0}\in \{90,100\}$ for $r_{0}$ ranging from 0% and 8%. We observe that the spreads between true and approximated prices depend on both parameters. However, it is noteworthy that the error remains under 1% in most cases.

Figure 2 Comparison of exact and PINN prices of a $T^{*}=T=8$ years GMAB with: $S_{0}=100$ , $\boldsymbol{\pi }=(0.50,0.25,0.25)$ , $(r_{0},\mu _{0},\lambda _{0})$ of Table 2. Upper plots: values at issuance and expiry, for $A_{t}\in [23.82,414.15]$ . Lower plots: values at issuance with $A_{0}=100$ and 90, for $r_{0}\in [0.00,0.08]$ .

To understand the circumstances in which relative errors are the highest, we analyze them by subgroups of contracts. In order to have a sufficient number of representatives in each sub-category, we generate a third validation set, this time consisting of 500,000 contracts instead of 100,000. Table 6 presents the averages and standard deviations of relative errors for 20 classes of contracts grouped by asset values. We observe that the highest errors occur with the lowest asset values. The left heatmap in Fig. 3 displays relative errors by categories of asset values and residual contract maturities, $T^{*}-t$ . We draw two conclusions from this graph. First, for $A_{t}\leq 80$ , the error increases with the residual contract duration and may be significant. Second, the errors in the usual pricing range, around the guarantee level $G_{m}=100$ , remain close to 0.01% regardless of the time horizon.

Table 6. Mean and standard deviations (std) of relative errors for contracts grouped by asset values

It is worth noting that for low asset values, GMAB is deeply out of the money, and the contract price is very close to that of a pure endowment for which we have a closed-form expression. One way to improve accuracy would be to add a third penalty term to the loss function (31). This term would penalize the quadratic spread between PINN and endowment prices for low asset values. However, we deliberately avoided this solution to observe how the network behaves when provided with minimal information.

Table 7. Mean and standard deviations (std) of relative errors for contracts grouped by interest rates

Figure 3 Left plot: heatmap of relative errors per asset and maturity classes. Right plot: heatmap of relative errors per interest rate and maturity classes.

Table 7 presents the means and standard deviations of errors for 10 groups of contracts, categorized by ranges of interest rates, $r_{t}$ . Regardless of the category, the average relative error remains close to 0.14%. The right heatmap in Fig. 3 suggests that, on average, errors increase with maturity. However, this observation should be qualified by the noticeable variability of errors between consecutive interest rate classes for the same maturity. This variability indicates that the increase in errors is partly attributed to other factors, such as the low asset values present within each of the interest rate classes.

Tables 8 and 9 provide information about errors for contracts categorized by classes of mortality rates. Similar to interest rates, the average error remains small and close to 0.14% across all class. The heatmaps of Fig. 4 show a slight increase in errors with the contract time to expiry. Once again, we observe a higher variability of errors for longer maturities. As mentioned earlier, this variability signals that the increase in errors is influenced by factors beyond mortality rates.

Table 8. Mean and standard deviations (std) of relative errors for contracts grouped by mortality rates $\mu _{x+t}$

Table 9. Mean and standard deviations (std) of relative errors for contracts grouped by hedging mortality rates $\lambda _{x+t}$

Figure 4 Left plot: heatmap of relative errors per $\mu _{x+t}$ (in $^{\circ }$ /oo) and maturity classes. Right plot: heatmap of relative errors per $\lambda _{x+t}$ (in $^{\circ }$ /oo) and maturity classes.

Tables 10 and 11 provide statistics on relative errors for contracts categorized by fractions of stocks and bonds (including interest rate and mortality bonds). The left plot in Fig. 5 displays a heatmap of errors based on these two parameters. Interestingly, we do not observe any particular trend, and the pricing error remains acceptable regardless of the asset mix. This indicates that the PINN can be used for ALM purposes, particularly in approximating the impact of changes in asset mix on contract fair values.

Table 10. Mean and standard deviations (std) of relative errors for contracts grouped by fraction of stocks

Table 11. Mean and standard deviations (std) of relative errors for contracts grouped by fraction of bonds

Figure 5 Heatmap of relative errors per $\pi _{S}$ and $(\pi _{P}+\pi _{D})$ asset allocation classes.

Nevertheless, the right heatmap of Fig. 5 reveals that the error increases with the time to expiry. To illustrate this, we compare the exact and PINN prices of a 5- and 8-year contract in Fig. 6, where $\pi _{S}$ ranges from 0 to 1 and $\pi _{D}=\pi _{M}=\frac{1-\pi _{S}}{2}$ . This graph shows that the PINN captures the trend of the relationship between prices and the fraction of stocks, but the error may reach 1% for certain combinations of parameters and risk factors. While this error may appear relatively high, it should be considered in the context of existing alternatives. For instance, in risk management, a common approach to valuing a contract at a future date relies on the least aquare Monte-Carlo (LSMC) method. As demonstrated in Hainaut & Akbaraly (Reference Hainaut and Akbaraly2023) for a similar contract, the average valuation error often exceeds 1% and can even surpass 3% in extreme cases. Replacing LSMC with the PINN in such a context would clearly reduce the overall valuation error.

Figure 6 Comparison of exact and PINN prices of a 5- and 8-year GMAB with: $S_{0}=100$ , $\pi _{S}\in [0,1]$ and $\pi _{D}=\pi _{M}=\frac{1-\pi _{S}}{2}$ . $(r_{0},\mu _{0},\lambda _{0})$ of Table 2. Upper plots: values at issuance and expiry, for $A_{t}\in [23.82,414.15]$ . Lower plots: values at issuance with $A_{0}=100$ and 90, for $r_{0}\in [0.00,0.08]$ .

We conclude from the error analysis that the average accuracy of a PINN is sufficient for risk management and ALM purposes. The significant advantage of using this approach is the substantial time savings once the calibration is completed. Computing prices for a sample set of 500,000 contracts, each with varying durations, asset mixes, and market conditions, takes less than 45 seconds. In contrast, evaluating the same portfolio using a standard MC method would require a considerable amount of time because we would need to regenerate a sufficient number of sample paths for each initial combination of risk factors. The PINN model is capable of pricing an impressive range of contracts, and with additional computing power, it is likely possible to further improve both accuracy and the diversity of contracts.