2.1 Financial market

Denote $(\Omega, \mathcal{F}, \mathbb{F}, {\mathbb P})$ as a complete filtered probability space, where ${\mathbb P}$ is a real-world probability measure and $\mathbb{F} \;:\!=\; \{\mathcal{F}_t\}_{t \geq 0}$ is a right-continuous and ${\mathbb P}$ -complete filtration, representing all the information up to time t. The stochastic processes in this paper are supposed to be well defined in the probability space and the moments of random variables are made under the probability measure ${\mathbb P}$ . In addition, we denote by ${\mathbb E} [\!\cdot\!]$ the expectation of the random variables under probability ${\mathbb P}$ .

The financial market consists of two assets, one risk-free asset $S_0(t)$ and one risky asset S(t). We assume that $S_0(t)$ is growing at a constant interest rate r, that is,

\begin{align*} \frac{{\mathop{}\!\mathrm{d}} S_0(t)}{S_0(t)} &= r {\mathop{}\!\mathrm{d}} t, \end{align*}

and S(t) is modeled by a general jump-diffusion model that allows for stochastic volatility and returns jumps (e.g., Bates, Reference Bates2000; Pan, Reference Pan2002 and Bégin, Reference Bégin2020),

(2.1) \begin{align} \frac{{\mathop{}\!\mathrm{d}} S(t)}{S(t-\!)} &= [r+ \lambda \times \nu(t) + \theta \times \nu(t) \times (\mu - \mu^{\mathbb{Q}})]{\mathop{}\!\mathrm{d}}{t} + \sqrt{\nu(t)}{\mathop{}\!\mathrm{d}} W_S(t) + {\mathop{}\!\mathrm{d}} \left( \sum_{n=1}^{N_t} (e^{Z_{n}} - 1) \right) - \mu \times \theta \times \nu(t) {\mathop{}\!\mathrm{d}} t, \nonumber \\[5pt] {\mathop{}\!\mathrm{d}} \nu(t) &= \kappa_{\nu}(\bar{\nu} - \nu(t)){\mathop{}\!\mathrm{d}} t + \sigma_{\nu} \sqrt{\nu(t)}\left[\rho_{\nu}{\mathop{}\!\mathrm{d}} W_S(t)+\sqrt{1-\rho_{\nu}^{2}}{\mathop{}\!\mathrm{d}} W_{\nu}(t)\right]. \end{align}

Here, $(W_S(t), W_{\nu}(t), W_l(t))$ is a three-dimensional standard Brownian Motion and we assume that $W_S(t), W_{\nu}(t)$ , and $W_l(t)$ are mutually independent. The instantaneous variance is modeled by $\nu(t)$ which is assumed to be a mean-reverting process with a long-run mean $\bar{\nu}$ , a mean-reverting rate of $\kappa_{\nu}$ , and a volatility coefficient of $\sigma_{\nu}$ , where $|\rho_{\nu}|\leq 1$ is a correlation coefficient between the volatility and the risky asset. We assume the parameters satisfy the condition $2\kappa_{\nu} \bar{\nu} \geq \sigma_{\nu}^2$ that ensures the positivity of $\nu(t)$ (see also Cox et al., Reference Cox, Ingersoll and Ross1985; Chen et al., Reference Chen, Nguyen and Stadje2018, and Chen et al., Reference Chen, Nguyen and Rach2021a for a similar condition). In addition, following Bégin (Reference Bégin2020), we include a jump component that follows a Poisson process $\{N_t\}_{t\geq 0}$ with a jump intensity of $\theta \nu(t)$ . The jump size process $\{Z_n\}_{n\in \mathbb{N}}$ follows a normal distribution with mean $\mu_z$ and standard deviation $\sigma_z$ , and we define $\mu = \mathbb{E}(\!\exp\!(Z_n)-1)$ as the expected jump return. The definition of the expected jump return under the risk-neutral measure $\mu^{\mathbb{Q}}$ is given in Section 2.2. Furthermore, we assume that the Poisson process $\{N_t\}_{t\geq 0}$ and the jump size process $\{Z_n\}_{n\in \mathbb{N}}$ are both independent of the Brownian motions.

The risk premium now consists of two components, $\lambda \nu(t)$ and $\theta \times \nu(t)\times (\mu - \mu^{\mathbb{Q} }) $ , corresponding to the premium for the return risks and the jump risks, respectively, and both are expressed as a constant multiplier of the volatility term $\nu(t)$ . We refer the readers to Pan (Reference Pan2002) for a thorough discussion on the specification of the risk premium and model calibration.

Furthermore, we model the salary uncertainty L(t) using an OU process, similar to Wang (Reference Wang2006; Reference Wang2009) and Luo (Reference Luo2017), that is,

(2.2) \begin{align} {\mathop{}\!\mathrm{d}} L(t)= \kappa_l(\bar{L}(t)-L(t)){\mathop{}\!\mathrm{d}} t + \sigma_l\sqrt{\nu(t)}\left[\rho_{lS}{\mathop{}\!\mathrm{d}} W_S(t)+\rho_{l\nu}{\mathop{}\!\mathrm{d}} W_{\nu}(t)+\sqrt{1-\rho_{lS}^{2}-\rho_{l\nu}^{2}}{\mathop{}\!\mathrm{d}} W_l(t)\right]. \end{align}

Here, $\kappa_l>0$ governs the speed of convergence to a deterministic function $\bar{L}(t)$ and we use an exponential function $\bar{L}(t) = \exp\!(\psi t)$ to roughly align with the average nominal salary growth over the past few decades.Footnote ⁵ $\sigma_l \times \sqrt{v(t)}$ is the volatility of the salary, and $\rho_{lS}$ and $\rho_{l\nu}$ are the correlation coefficients between the salary and the risky asset and between the salary and the volatility, respectively. Without loss of generality, we unitize the salary such that $L(0) = 1$ . Alternatively, L(t) can be regarded as the inflation index, which reflects the average cost of living of retirees. In this paper, we do not differentiate the salary and the inflation risks.

2.2 Risk-neutral dynamics and volatility index

The estimate of the risk premium and the benchmark $\overline{\text{VIX}^2}$ involves the specification of the risk-neutral dynamics. Due to the additional sources of uncertainty (i.e., the stochastic volatility and the random jump sizes), the market is incomplete in our setting. As a result, the state-price density is not unique. We follow the method from Pan (Reference Pan2002) by focusing on a candidate state-price density that prices the two important sources of risks: diffusive price shocks and jump risks. Our candidate state-price density is specified in Appendix B.1. Under the risk-neutral measure $\mathbb{Q}$ associated with the candidate state-price density, the dynamics of the risky asset are in the following.

(2.3) \begin{align} \frac{{\mathop{}\!\mathrm{d}} S(t)}{S(t-\!)} &= r {\mathop{}\!\mathrm{d}} t + \sqrt{\nu(t)} {\mathop{}\!\mathrm{d}} W_S^{\mathbb{Q} }(t) + {\mathop{}\!\mathrm{d}} \left(\sum_{n=1}^{N_t^{\mathbb{Q} } } \left(e^{Z_n^{\mathbb{Q} }}-1 \right) \right) - \mu^{\mathbb{Q} } \theta \nu(t) {\mathop{}\!\mathrm{d}} t\nonumber \\[5pt] {\mathop{}\!\mathrm{d}} \nu(t) &= \kappa_{\nu} (\bar{\nu} - \nu(t) ){\mathop{}\!\mathrm{d}} t + \sigma_{\nu} \sqrt{\nu(t)} \left[ \rho_{\nu} {\mathop{}\!\mathrm{d}} W_S^{\mathbb{Q}}(t) + \sqrt{1-\rho_{\nu}^2} {\mathop{}\!\mathrm{d}} W_{\nu}^{\mathbb{Q} }(t) \right], \end{align}

where $\mu^{\mathbb{Q} } = \mathbb{E}^{\mathbb{Q} } \left( \!\exp\!(Z^{\mathbb{Q}}_n) - 1 \right)$ . Here, $\left(W^{\mathbb{Q}}_S(t), W^{\mathbb{Q}}_{\nu}(t), W^{\mathbb{Q}}_l(t)\right)$ , $\{N^{\mathbb{Q}}_t\}_{t\geq 0}$ , and $\{Z^{\mathbb{Q}}_n\}_{n\in \mathbb{N}}$ are a three-dimensional standard Brownian Motion, the Poisson process with jump intensity of $\theta\cdot \nu(t)$ , and the jump size process with mean $\mu_z^{\mathbb{Q}}$ and variance $\sigma_z^2$ under $\mathbb{Q}$ , respectively.Footnote ⁶

Given the risk-neutral dynamic, the volatility index that measures the market’s expectation of the annualized volatility over the next 30 days has been derived explicitly by Lin (Reference Lin2007). Here, we present the results only and refer interested readers to Lin (Reference Lin2007) for details. The squared volatility index $\text{VIX}_t^2$ can be expressed as

\begin{align*} \text{VIX}_t^2 &= a_{\text{VIX}} \times \nu(t) + b_{\text{VIX}}, \end{align*}

where

\begin{align*} a_{\text{VIX}} &= 100^2\times \left( 1 + 2\theta\cdot (\mu^{\mathbb{Q} } - \mu_z^{\mathbb{Q} }) \right) \cdot \frac{1 -e^{{\kappa_{\nu}} \cdot \tau }}{{\kappa_{\nu}} \cdot \tau },\\[5pt] b_{\text{VIX}} &=100^2\times \left( 1 + 2\theta\cdot (\mu^{\mathbb{Q} } - \mu_z^{\mathbb{Q} }) \right)\cdot {\bar{\nu}} \cdot \left[ 1-\frac{1 -e^{{\kappa_{\nu}} \cdot \tau }}{{\kappa_{\nu}} \cdot \tau } \right], \end{align*}

and $\tau = \frac{30}{365}$ .

The benchmark used for volatility adjustment in Bégin (Reference Bégin2020) is the asymptotic mean of $\text{VIX}^2_t$ under the risk-neutral measure, that is,

\begin{align*} \overline{\text{VIX}^2} = \lim_{t\rightarrow \infty} \mathbb{E}^{\mathbb{Q}}[\text{VIX}_t^2] = a_{\text{VIX}} \times \bar{\nu} + b_{\text{VIX}}. \end{align*}

2.3 Population and pension plan provision

In this paper, we consider a continuous population where all employees join the work at age A and retire at age R with a maximum attainable age of $\omega$ . We assume there is no other pre-exit of the pension plan except death. Denote $n_y(t)$ as the number of members aged y at time t, which is assumed to be a deterministic function of time t. We use the standard actuarial notation ${}_tp_y$ to represent the probability for a person aged y that survives for the next t years. We do not consider the longevity risk and assume the mortality risk is fully diversifiable, that is, $n_y(t) = n_A(t-(y-A)) \times {}_{y-A}p_A$ .

The pension plan is assumed to be fully funded by the existing members, where a fixed contribution rate of c applies to all active members. The benefit is quoted in terms of the instantaneous replacement rate b(t), which is not a guaranteed amount.Footnote ⁷ Later in this paper, we will include a discussion on the contribution adjustment (CA) scheme and the hybrid scheme that have been studied in Cui et al. (Reference Cui, De Jong and Ponds2011) and Bégin (Reference Bégin2020). The aggregate contribution and the aggregate benefit at time t are defined as

\begin{align*} &\text{Aggregate contribution}(t) = \int_A^{R} n_y(t)\times c \times L(t) \times e^{\eta_a (y-A)} {\mathop{}\!\mathrm{d}}{y} = c\times L(t) \times \mathcal{A}(t),\\[5pt] &\text{Aggregate benefit}(t) = \int_A^{\omega} n_y(t)\times b(t) \times L(t) \times e^{\eta_a(R-A)} {\mathop{}\!\mathrm{d}}{y} = b(t) \times L(t) \times \mathcal{R}(t), \end{align*}

where $\eta_a\geq 0$ is the promotional adjustment such that employees with longer career life generally have a higher salary. When $\eta_a = 0$ , the aggregate contribution and benefits are simply multipliers of the population sizes of the active workers $ \mathcal{A}(t)$ and the retirees $\mathcal{R}(t)$ , respectively.

2.4 Pension asset and liability

Given the financial models and the aggregate contributions and benefits, the pension asset X(t) has the following dynamics:

(2.4) \begin{align} {\mathop{}\!\mathrm{d}} X(t) &= (X(t)-\pi(t))\frac{{\mathop{}\!\mathrm{d}} S_0(t)}{S_0(t)} + \pi(t)\frac{{\mathop{}\!\mathrm{d}} S(t)}{S(t-\!)} + c\times L(t) \times \mathcal{A}(t){\mathop{}\!\mathrm{d}} t - b(t) \times L(t) \times \mathcal{R}(t){\mathop{}\!\mathrm{d}} t \nonumber\\[5pt] &= \left[rX(t)+\pi(t)\nu(t)(\lambda - \theta\mu^{\mathbb{Q} } )+(c\mathcal{A}(t)-b(t)\mathcal{R}(t))\times L(t)\right] {\mathop{}\!\mathrm{d}}{t} \nonumber\\[5pt] &\qquad +\pi(t)\sqrt{\nu(t)}{\mathop{}\!\mathrm{d}} W_{S}(t)+ \pi(t){\mathop{}\!\mathrm{d}} \left( \sum_{n=1}^{N_t} (e^{Z_{n}} - 1) \right), \end{align}

with the initial asset level $X(0) = x_0$ . The amount invested in the risky asset $\pi(t)$ and the replacement rate b(t) are control variables that are dynamically decided by the pension sponsor.

In this paper, we use the traditional unit credit (TUC) method to valuate the actuarial liability (AL). Since the benefit level is not guaranteed, we apply the TUC to a benchmark benefit level $\hat{b}$ such that

\begin{align*} \text{AL}(t) &= \int_{A}^{\omega} n_y(t) \times \underset{ \text{= Benefit Accrual}}{\underbrace{\frac{\min\!(R,y)-A}{R-A}\times \hat{b}}}\times L(t)\times e^{\eta_a(\!\min\!(R,y)-A)}\times \int_{\max (0, R-y)}^{\omega-y} {}_{s}p_{y} \times e^{-\phi s}{\mathop{}\!\mathrm{d}}{s} {\mathop{}\!\mathrm{d}}{y}\\[5pt] &= \mathcal{H}(t) \times L(t), \end{align*}

where $\phi$ is the actuarial discount rate and $\mathcal{H}(t)$ represents the AL value at time t in the real term.

2.5 Objective

We search for a welfare-enhancing risk-sharing design. We follow the standard approach to define the welfare function as the aggregate discounted expected utility function for all members during a planning horizon of T years. To avoid the situation where the interests of the future generations beyond the time T are sacrificed for the benefits of the current and nearby generations, we include a penalty term that depends on the terminal asset level X(T).

The objective is to maximize the welfare function by controlling the investment and the benefit payment strategies. Define the value function $H(t, x, \nu, l)$ as

(2.5) \begin{align} H(t,x,\nu,l) = \sup_{(\pi, b)\in\Pi} \mathbb{E}_{t,x,\nu,l}\left[\int_{t}^{T} e^{-\zeta (s - t)}\times \mathcal{R}(s) \times U( b(s)\times L(s);\; \gamma_r) {\mathop{}\!\mathrm{d}} s + \varrho \times e^{-\zeta (T-t)} \times U(X(T);\; \gamma_T) \right], \end{align}

where $E_{t,x,\nu,l} \;:\!=\;E[ \cdot |t, X(t) = x,\nu(t) = \nu,L(t) = l]$ , $\zeta$ is the discount rate for the time preference, $\varrho$ is the preference weight given to the utility function of the terminal asset level, and $\Pi$ is the admissible set defined in Appendix B.2. We adopt an exponential utility function such that $U(x;\;\gamma) = - \frac{1}{\gamma} \exp\!(\!-\!\gamma x)$ , where $\gamma$ is the risk aversion parameter. Note that we specify different values of the risk aversion parameters for the individual retirement benefit and the asset level at the terminal time, namely $\gamma_r$ and $\gamma_T$ .

3.1. Optimal benefit payment structure

To solve Problem (2.5), we use the Hamiltonian–Jacobi–Bellman approach and derive the explicit solution. The explicit solution and the optimal benefit payment strategy are summarized in the following theorem.

Theorem 3.1. The optimal instantaneous replacement rate $b^*(t)$ is

(3.1) \begin{align} b^{*}(t)=-{\displaystyle\frac{1}{\gamma_r L(t)}}\ln{\displaystyle\frac{\gamma_T\varrho A(t)}{\gamma_r}}+{\displaystyle\frac{\gamma_T}{\gamma_r L(t) }}[A(t)X(t)+\bar{A}(t)\nu(t)+\hat{A}(t)L(t)+\tilde{A}(t)], \end{align}

the optimal investment strategy $\pi^*(t)$ satisfies the following equation:

(3.2) \begin{align} \pi^{*}(t) \cdot \gamma_{T}A(t)-\theta\mathbb{E}[{\rm e}^{-\gamma_{T}A(t)\pi^{*}(t)\cdot({\rm e}^{Z_{n}}-1)}({\rm e}^{Z_{n}}-1)] = (\lambda-\theta\mu^{\mathbb{Q}})-\sigma_{\nu}\rho_{\nu}\gamma_{T}\bar{A}(t)-\sigma_l\rho_{lS}\gamma_{T}\hat{A}(t), \end{align}

and the corresponding value function $H(t,x,\nu,l)$ has an exponential form such that

\[ H(t,x,\nu,l)=-\frac{\varrho}{\gamma_T} {\rm e}^{-\gamma_{T}[A(t)x+\bar{A}(t)\nu+\hat{A}(t)l+\tilde{A}(t)]}. \]

Here A(t), $\bar{A}(t)$ , $\hat{A}(t),$ and $\tilde{A}(t)$ are the solution of the following system of differential equations:

(3.3) \begin{align} \left\{ \begin{array}{l} A_{t}(t) = -rA(t)+{\displaystyle\frac{\mathcal{R}(t)\gamma_{T}A^{2}(t)}{\gamma_{r}}}\\[5pt] \hat{A}_{t}(t)=-c\mathcal{A}(t)A(t)+\kappa_l\hat{A}(t)+{\displaystyle\frac{\mathcal{R}(t)\gamma_{T}A(t)\hat{A}(t)}{\gamma_{r}}}\\[5pt] \bar{A}_{t}(t)=-\pi^{*}(t)(\lambda-\theta\mu^{\mathbb{Q}})A(t)+\kappa_{\nu}\bar{A}(t)+{\displaystyle\frac{1}{2}}(\pi^{*})^{2}\gamma_{T}A^{2}(t)+{\displaystyle\frac{1}{2}}\sigma_{\nu}^{2}\gamma_{T}\bar{A}^{2}(t)+{\displaystyle\frac{1}{2}}\sigma_l^{2}\gamma_{T}\hat{A}^{2}(t)\\[5pt] \qquad \qquad +\pi^{*}(t)\sigma_{\nu}\rho_{\nu}\gamma_{T}A(t)\bar{A}(t)+\pi^{*}(t)\sigma_l\rho_{lS}\gamma_{T}A(t)\hat{A}(t) +\sigma_l\sigma_{\nu}(\rho_{lS}\rho_{\nu}+\rho_{l\nu}\sqrt{1-\rho_{\nu}^{2}})\gamma_{T}\bar{A}(t)\hat{A}(t)\\[5pt] \qquad \qquad +{\displaystyle\frac{\theta}{\gamma_{T}}}\mathbb{E}[{\rm e}^{-\gamma_{T}A(t)\pi^{*}(t)\cdot ({\rm e}^{Z_{n}}-1)}-1] +{\displaystyle\frac{\mathcal{R}(t)\gamma_{T}A(t)\bar{A}(t)}{\gamma_{r}}}\\[5pt] \tilde{A}_{t}(t)=-{\displaystyle\frac{\zeta}{\gamma_{T}}}-\kappa_{\nu}\bar{\nu}\bar{A}(t)-\kappa_l\bar{L}(t)\varrho\hat{A}(t)-{\displaystyle\frac{\mathcal{R}(t)A(t)}{\gamma_{r}}}\ln(\varrho A(t))+{\displaystyle\frac{\mathcal{R}(t)\gamma_{T}A(t)\tilde{A}(t)}{\gamma_{r}}}+{\displaystyle\frac{\mathcal{R}(t)A(t)}{\gamma_{r}}}, \end{array}\right. \end{align}

with terminal conditions $A(T)=1$ and $\bar{A}(T)=\hat{A}(T)=\tilde{A}(T)=0$ .Footnote ⁸

The solution for the optimal investment strategy aligns well with other stochastic control studies using an exponential utility function (e.g., Merton, Reference Merton1969; Yang and Zhang, Reference Yang and Zhang2005, and Wang, Reference Wang2007). More importantly, we focus on the risk-sharing structure of the benefit payment. Rearranging Equation (3.1), we have the following important remark.

Remark 1. The optimal benefit payment can be written in the following structure:

(3.4) \begin{align} \text{benefit}(t) = b^*(t)\cdot L(t)=\bar{b}(t) + \beta_A(t) \frac{X(t) - \xi_A(t) \times \text{AL}(t) }{\mathcal{R}(t)} - \beta_{\text{VIX}}(t) \frac{\text{VIX}^2_t - \xi_{\text{VIX}}(t)\times \overline{\text{VIX}^2}}{\mathcal{R}(t)}, \end{align}

where the risk-sharing parameters $\beta_A(t)$ and $\beta_{\text{VIX}}$ are defined as

\begin{align*} \beta_A(t) &= \frac{\gamma_T \mathcal{R}(t) A(t)}{\gamma_r},\qquad \beta_{\text{VIX}}(t) = \frac{-\gamma_T \mathcal{R}(t) \bar{A}(t)}{\gamma_r a_{\text{VIX}} }, \end{align*}

and the values for $\xi_A(t)$ and $\xi_{\text{VIX}}(t)$ can be arbitrarily chosen such that the TB level $\bar{b}(t)$ is defined as

\begin{align*} \bar{b}(t) &= -\frac{1}{\gamma_r}\ln \left(\frac{\gamma_T}{\gamma_r}\varrho A(t)\right) + \frac{\gamma_T}{\gamma_r} \tilde{A}(t) - \frac{\gamma_T \bar{A}(t)}{\gamma_r a_{\text{VIX}} } (b_{\text{VIX}}-\xi_{\text{VIX}}(t)\times \overline{\text{VIX}^2} )+ \frac{\gamma_T}{\gamma_r} \left( A(t) \xi_A(t) \mathcal{H}(t) + \hat{A}(t) \right) L(t). \end{align*}

In general, the structure of the optimal benefit payment (a.k.a. Equation (3.4)) coincides with the proposal made by Bégin (Reference Bégin2020), except that all the risk-sharing parameters (i.e., $\beta_A(t)$ and $\beta_{\text{VIX}}(t)$ ) in Equation (3.4) are deterministic functions of time t instead of constants when the planning time horizon is finite or when the population structure changes in the future. Therefore, the implication of the optimal structure is straightforward, that is, the actual benefit ( $b(t) \cdot L(t)$ ) deviates from the target level $\bar{b}(t)$ when the asset level is different from a certain threshold ( $\xi_A(t)$ ) of the AL or when the VIX index is different from a certain threshold ( $\xi_{\text{VIX}}$ ) of the benchmark $\overline{\text{VIX}^2}$ . To further align with the proposal in Bégin (Reference Bégin2020) and simplify the following numerical analysis, we set $\xi_A(t) = \xi_{\text{VIX}} = 1$ to represent that the risk-sharing threshold is 100% of the AL and the volatility benchmark. We keep the terms $\xi_A(t)$ and $\xi_{\text{VIX}}(t)$ in Equation (3.4) to stress that these thresholds are allowed to be different from 100%.Footnote ⁹ In what follows, we examine each of the three components in Equation (3.4), namely the performance adjustment $\beta_A$ , the volatility adjustment $\beta_{\text{VIX}}$ , and the TB level $\bar{b}$ .

3.2. Performance adjustment $\beta_A$

The performance adjustment $\beta_A(t)$ can be explicitly expressed as

(3.5) \begin{align} \beta_A(t) = \frac{ \mathcal{R}(t) } { \frac{\gamma_r}{\gamma_T} e^{-(T - t) r}+ \int_t^T e^{-(s-t)r} \mathcal{R}(s) {\mathop{}\!\mathrm{d}}{s}}, \end{align}

which is independent of most of the economic parameters and guaranteed to be positive. In fact, it is only affected by the choice of the risk-free rate r, the ratio between the risk aversion parameters $\frac{\gamma_T}{\gamma_r}$ , and the size of the retired population $\mathcal{R}(t)$ . The effect of the risk aversion ratio $\frac{\gamma_T}{\gamma_r}$ becomes insignificant for a sufficiently long horizon. In particular, we have

(3.6) \begin{align} \lim_{T\rightarrow \infty}\beta_A(t;T) =\frac{ \mathcal{R}(t) }{ \int_{t}^{\infty} e^{-(s-t)r} \mathcal{R}(s) {\mathop{}\!\mathrm{d}}{s} }. \end{align}

We can observe from Equation (3.6) that $\beta_A$ only depends on the relative population size of the current retirees compared with that of the future retirees. If we have a fast-growing population, then future generations have more ability in absorbing the pension deficit. Therefore, it is optimal to stabilize the retirement income for current retirees with more deficit transferred to the future, which results in a smaller $\beta_A$ . Notice that the discount rate for the time preference (or intergenerational discount rate) $\zeta$ is not involved in Equation (3.6), which means that the sponsor does not need to worry about the issues of intergenerational fairness when determining the performance adjustment.

It is important to emphasize that the performance adjustment $\beta_A$ is determined independently of the stochastic volatility assumption. We can show that $\beta_A(t)$ still has the same formula when the risky asset is assumed to follow a geometric Brownian motion as Merton (Reference Merton1969).Footnote ¹⁰ This suggests that the performance adjustment term is model-robust. The volatility adjustment can be regarded as a supplementary adjustment on top of the performance adjustment.

3.3. Volatility adjustment $\beta_{\text{VIX}}$

We first show that the volatility adjustment parameter $\beta_{\text{VIX}}(t)$ is guaranteed to be non-positive when the jump and salary risks are excluded. This result is summarized in the following remark.

Remark 2. If the jump and salary risks are not considered, then

(3.7) \begin{equation} \beta_{\text{VIX}}(t)=\left\{\begin{array}{ll} {\displaystyle -\frac{\gamma_T}{\gamma_r}\frac{\mathcal{R}(t)}{a_{\text{VIX}}} \frac{\nu_{1}-\nu_{1}{\rm e}^{-\frac{\sigma_{\nu}^{2}(1-\rho_{\nu}^{2})\gamma_{T}(\nu_{1}-\nu_{2})}{2}(T-t)}}{1-\frac{\nu_{1}}{\nu_{2}}{\rm e}^{-\frac{\sigma_{\nu}^{2}(1-\rho_{\nu}^{2})\gamma_{T}(\nu_{1}-\nu_{2})}{2}(T-t)}}}, & \rho_{\nu}\neq\pm1,\\[5pt] {\displaystyle -\frac{\mathcal{R}(t)}{a_{\text{VIX}}} \frac{\lambda^{2}}{2\gamma_{r}}}\int_{t}^{T}{\rm e}^{-\int_{0}^{w}(\lambda\sigma_{\nu}+\kappa_{\nu}+\beta_{A}(s))\cdot(1+\mathbb{1}_{s\leq t}) {\rm d}s}{\rm d}w, & \rho_{\nu}=1,\;\;\lambda\sigma_{\nu}\rho_{\nu}+ \kappa_{\nu}\neq-\beta_{A}(t), \\[5pt] {\displaystyle -\frac{\mathcal{R}(t)}{a_{\text{VIX}}} \frac{\lambda^{2}}{2\gamma_{r}}} \int_{t}^{T}{\rm e}^{-\int_{0}^{w}(-\lambda\sigma_{\nu}+\kappa_{\nu}+\beta_{A}(s))\cdot(1+\mathbb{1}_{s\leq t}){\rm d}s}{\rm d}w, & \rho_{\nu}=-1,\;\;\lambda\sigma_{\nu}\rho_{\nu}+ \kappa_{\nu}\neq-\beta_{A}(t), \\[5pt] {\displaystyle -\frac{\mathcal{R}(t)}{a_{\text{VIX}}} \frac{\lambda^{2}}{2\gamma_{r}}} (T-t), & \rho_{\nu}=\pm1,\;\;\lambda\sigma_{\nu}\rho_{\nu}+ \kappa_{\nu}=-\beta_{A}(t), \end{array} \right. \end{equation}

where

\begin{equation*} \nu_{1,2}={\displaystyle\frac{-\lambda\sigma_{\nu}\rho_{\nu}- \kappa_{\nu}-\beta_{A}(t)\pm\sqrt{(\lambda\sigma_{\nu}\rho_{\nu}+ \kappa_{\nu}+\beta_{A}(t))^{2}+\lambda^{2}\sigma_{\nu}^{2}(1-\rho_{\nu}^{2})}} {\sigma_{\nu}^{2}(1-\rho_{\nu}^{2})\gamma_{T}}}. \end{equation*}

It is clear from Remark 2 that

\begin{align*} \beta_{\text{VIX}}(t) \leq 0 \end{align*}

by observing that $\nu_{1} \geq 0\geq \nu_2$ . The derivation of Equation (3.7) can be found in Appendix C.2. We can also show that $\beta_{\text{VIX}}(t)$ is non-positive numerically when both the jump and the salary risks are included in Section 4 (see Figures 2 and 3).

The non-positivity of $\beta_{\text{VIX}}(t)$ indicates that the actual benefit will be increased when the volatility index is higher than the long-term benchmark (i.e., $\text{VIX}^2>\overline{\text{VIX}^2}$ ) if all other items are kept unchanged in Equation (3.4). It can be explained by the positive choices in the risk premium $\lambda$ , where a higher risk premium increases the equity’s investment return and raises the benefit level further. Interestingly, the findings imply that retirees will receive compensation for potential benefit reductions during a market crash where the value of the pension asset is likely to decline sharply and the volatility index is typically quite high. This observation suggests a counter-cyclical risk-sharing policy and is in line with Chen et al. (Reference Chen, Li, Wang and Zhu2023). Notice that the non-positive values of $\beta_{\text{VIX}}(t)$ are the direct result of our optimization problem. This is in contrast with the non-negative assumption made in Bégin (Reference Bégin2020). While Bégin (Reference Bégin2020) sets the assumption due to financial fairness, our welfare-maximizing design provides a different interpretation on the volatility adjustment $\beta_{\text{VIX}}(t)$ .

3.4. Partial indexation of TB $\bar{b}$

The TB level $\bar{b}(t)$ in Remark 1 can be decomposed in the following:

\begin{align*} \bar{b}(t) = \bar{b}_F(t) + \bar{b}_L(t) \times L(t), \end{align*}

where

\begin{equation*} \bar{b}_F(t) = -\frac{1}{\gamma_r}\ln \left(\frac{\gamma_T}{\gamma_r}\varrho A(t)\right) + \frac{\gamma_T}{\gamma_r} \tilde{A}(t) - \frac{\gamma_T \bar{A}(t)}{\gamma_r a_{\text{VIX}} } \left(b_{\text{VIX}}-\xi_{\text{VIX}}(t)\times \overline{\text{VIX}^2} \right) \end{equation*}

represents the fixed amount of the TB that is independent of the salary risk, and

\begin{equation*} \bar{b}_L(t) = \frac{\gamma_T}{\gamma_r} \left( A(t) \xi_A(t) \mathcal{H}(t) + \hat{A}(t) \right) \end{equation*}

stands for the variable amount that reflects the COLA. The benefit is fully indexed to a wage index when $\bar{b}_F = 0$ , and no indexation is applied when $\bar{b}_L = 0$ . It is evident that both $\bar{b}_F$ and $\bar{b}_L$ are strictly positive when $\xi_{A} = \xi_{\text{VIX}} = 1$ . This indicates that the optimal retirement benefit adopts a partial COLA, which aligns well with the market practice for many existing DB plans.Footnote ¹¹ In contrast, most of the literature finds either full indexation or no indexation.

3.5. Hybrid scheme with adjustable contribution rate c(t)

In this section, we include a short discussion on the optimal CA. Specifically, we now set the contribution rate as an additional control variable c(t) and the objective is still to maximize the welfare function. Active members’ utility is based on their disposable income after the contribution, and an exponential utility function with risk aversion parameter $\gamma_a$ is used. Here, we summarize the results in the following remark; the details of the problem formulation along with the derivation are deferred to Appendix D.

Remark 3. The optimal contribution $c(t)\cdot L(t)$ and benefit payment $b(t)\cdot L(t)$ can be written in a similar structure as in Bégin (Reference Bégin2020),

\begin{align*} \text{Contribution}(t) &= \bar{c}(t) - \alpha_A(t) \frac{X(t) - \text{AL}(t) }{\mathcal{A}(t) } + \alpha_{\text{VIX}}(t) \frac{\text{VIX}^2_t - \overline{\text{VIX}^2} }{\mathcal{A}(t) },\\[5pt] \text{Benefit}(t) &= \bar{b}(t) + \beta_A(t) \frac{X(t) - \text{AL}(t) }{\mathcal{R}(t) } - \beta_{\text{VIX}}(t) \frac{\text{VIX}^2_t - \overline{\text{VIX}^2} }{\mathcal{R}(t) }, \end{align*}

where the risk-sharing parameters $\alpha$ and $\beta$ and the target contribution $\bar{c}(t)$ and TB $\bar{b}(t)$ are defined in Appendix D with the same interpretation and the same structure as $\beta$ in Remark 2.

The remark justifies the linear risk-sharing proposals made by Cui et al. (Reference Cui, De Jong and Ponds2011) and Bégin (Reference Bégin2020) on the optimal hybrid design. Most observations on the TB plan remain valid for the hybrid plan and more importantly, the relationship of the risk-sharing parameters between the contribution and the benefit exhibits the following constant ratio:

\begin{align*} \frac{\alpha_A(t)}{\beta_A(t) }= \frac{\alpha_{\text{VIX}}(t) }{\beta_{\text{VIX}}(t) } = \frac{\gamma_r \mathcal{A}(t) }{ \gamma_a \mathcal{R}(t) }. \end{align*}

This indicates the values for $\alpha$ and $\beta$ differ only due to the difference in the population size and risk aversion between active workers and retirees. In particular, when $\gamma_r = \gamma_a$ , the amount of risk borne by any group is proportional to the group population size.Footnote ¹² In the usual circumstance where $\mathcal{A}(t)>\mathcal{R}(t)$ and $\gamma_a < \gamma_r$ , active workers will carry relative higher risks than retirees. Moreover, the constant ratio relationship implies that $\alpha_{\text{VIX}}$ should be non-positive, in contrast with the non-negative assumption used in Bégin (Reference Bégin2020). In the end, when $\gamma_a = \gamma_r$ it can be shown that

\begin{align*} \text{contribution}(t) + \text{benefit}(t) = L(t). \end{align*}

This means that the after-contribution income for the active worker is the same as the retirement income for the retirees, which implies a 100% replacement rate. Clearly, this is not the general consensus for a pension design as a replacement rate around 60% is often deemed to be reasonable for a lifetime worker and therefore it would be reasonable to set $\gamma_a\approx 60\% \times \gamma_r$ .

4.1. Performance adjustment $\beta_A$

Due to the decreasing relative risk aversion of the exponential utility function, it is appropriate to set different values for the risk aversion parameters of the retirement income and the terminal asset level because they are significantly different in scale.Footnote ¹⁵ However, the choices of $\gamma_r$ and $\gamma_T$ are not straightforward because (1) they represent two different quantities that are not comparable, with the former one focusing on an individual’s income adequacy while the latter one on the overall pension solvency and (2) both the retirement income and the terminal asset level are random. Therefore, we use an alternative approach to provide a new interpretation between $\gamma_r$ and $\gamma_T$ . Recall that $\beta_A(t)$ represents the percentage of deficit/surplus that will be distributed to the retirees and that $\beta_A(T) = \frac{\gamma_T \mathcal{R}(T) }{\gamma_r}$ . The choice of $\frac{\gamma_T}{\gamma_r}$ has a natural upper limit of $\frac{1}{\mathcal{R}(T)}$ so that the pension fund cannot distribute more than 100% of the deficit to the retirees. Although the choices of $\gamma_r$ and $\gamma_T$ under this approach deviate from their original definitions with respect to utility functions, they reconcile with the objective of balancing the conflicting interests between different generations. The value of $\beta_A$ (and hence the ratio of $\gamma_T$ over $\gamma_r$ ) represents the level of risk bored by the current generations and thus indirectly reflects the risk transfer to the future.

Figure 1 displays the effect of $\frac{\gamma_T}{\gamma_r}$ (expressed in term of $\beta_A(T)$ within the range 1–100%) on $\beta_A(t)$ . We select $\beta_A(T) = 3\%$ in our benchmark scenario to match the optimal value obtained in Bégin (Reference Bégin2020). Note that this value also provides the most time-stable performance adjustment as shown in the graph. This figure demonstrates that $\beta_A$ is relatively stable for all four cases (at least for the initial period of roughly twenty years). The shape of $\beta_A(t)$ is usually flat over the initial period of time and increases/decreases sharply when the time approaches T (with speed controlled by the risk-free rate and the population growth rate of the retirees).

Figure 1. $\beta_A(t)$ under different $\beta_A(T) = \frac{\gamma_T}{\gamma_r }\mathcal{R}(T) $ .

Figure 2. $\beta_{\text{VIX}}(t)$ under different $\beta_A(T) = \frac{\gamma_T}{\gamma_r} \mathcal{R}(T) $ .

Figure 3. $\beta_{\text{VIX}}(0)$ under different $\lambda$ .

4.2. Volatility adjustment $\beta_{\text{VIX}}$

Figure 2 demonstrates the values of the volatility adjustment across time under different values of $\gamma_T$ (expressed in terms of $\beta_A(T)$ ). First, we can observe that $\beta_{\text{VIX}}$ appears to be more stable over a longer period of time than $\beta_A(t)$ for all four cases. It stays flat for over 35 years before a sharp convergence to 0 near maturity. Second, the volatility adjustment is robust with respect to different choices of $\gamma_T$ within a reasonable range, which is similar to the performance adjustment. Third, the volatility adjustment $\beta_{\text{VIX}}$ is negative and significantly differs away from zero for all the cases, which echoes Remark 2 where the jump and salary risks are not considered. This indicates that the volatility adjustment in fact possesses the features of counter-cyclical risk-sharing policies.

Figure 3 displays the value of $\beta_{\text{VIX}}$ with different risk premia ( $\lambda$ ) and different correlation coefficients between the volatility and the risky asset ( $\rho_{\nu}$ ) at time $t=0$ . We can observe that the effect of volatility adjustment is most significant (i.e., largest in the absolute term $|\beta_{\text{VIX}}|$ ) when the risk premia are high and $\rho_{\nu}$ is highly negative, which highlights the hedging effect of the volatility adjustment term during the bear market. When we move from negative correlation coefficients to positive, we can observe that the hedging effect of the volatility adjustment diminishes since the high-risk premia will offset the market crash. In addition, the higher the risk premia $\lambda$ , the better the expected investment returns will be achieved, and thus more surplus are expected to be available for distribution to the retirees.

With a non-negative constraint, Bégin (Reference Bégin2020) obtain obtain an optimal benefit structure that does not have a volatility adjustment ( $\beta_{\text{VIX}}$ = 0).Footnote ¹⁶ We deem that this is partly because of misspecification in the risk premia. If we adopt the same assumption as Bégin (Reference Bégin2020) such that the risk premia are independent of the volatility term (i.e., $\lambda = 0$ ), then we would also obtain a zero volatility adjustment as shown in Figure 3. However, as demonstrated by Pan (Reference Pan2002), investors seek a higher expected return on riskier assets, and therefore a large positive $\lambda$ is more consistent with the market expectation. In such a case, the optimal risk-sharing design should always contain a negative $\beta_{\text{VIX}}$ in its benefit adjustment (BA).

Figure 4 illustrates the historical BA based on the optimal $\beta_A(t)$ and $\beta_{\text{VIX}}$ using the monthly S&P 500 data from 2006 to 2014. The equity market has done poorly during this time (i.e., behind inflation), and BAs have been persistently negative. All terms are expressed in real values (by dividing the average wage index) and we adopt a 50-50 investment strategy in risky and risk-free assets in this illustration. We can observe that the performance adjustment is dominating the overall adjustment and the volatility adjustment becomes significant only during the 2008 financial crisis. The volatility adjustment is clearly counter-cyclical to provide a hedge on top of the performance adjustment, which is pro-cyclical. Indeed, we find that the variation in the benefit adjustment (i.e., the standard deviation of $\text{Adjustment}_{t+\Delta} - \text{Adjustment}_t$ ) is reduced by 13% when the volatility adjustment is introduced and the adjustment volatility (i.e., the standard deviation of BA) is reduced by 9%.

Figure 4. Sample of historical benefit adjustments.

4.3. Partial indexation of the optimal benefit

To measure the partial indexation level of the TB, recall that given $L(0) = 1$ we have

\begin{align*} \bar{b}(t) = \bar{b}_F(t) + \bar{b}_L(t) \cdot L(t) = \left(\bar{b}_F(t) + \bar{b}_L(t)\right) \cdot \left[ 1 + \frac{\bar{b}_L(t) }{\bar{b}_F(t) + \bar{b}_L(t)} \cdot \left( L(t) - 1 \right) \right].\end{align*}

Since $L(t) - 1$ is the cumulative inflation growth up to time t, we define $\mathcal{I}(t) = \frac{\bar{b}_L(t) }{\bar{b}_F(t) + \bar{b}_L(t)}$ as the partial indexation rate.Footnote ¹⁷

Figure 5 plots the levels of partial indexation at time $t=0$ over different $\varrho$ (weight given to the terminal asset) and $\zeta$ (time preference). Since $\varrho$ and $\zeta$ only appear in the ordinary differential equation (ODE) of $\tilde{A}(t)$ and thus do not affect the risk-sharing parameters $\beta$ or $\bar{b}_L$ , Figure 5 also illustrates the sensitivity of $\bar{b}_F$ with respect to $\varrho$ and $\zeta$ . We first find that the indexation level is within the acceptable range for DB practice, that is, between 59% and 63%, and is comparatively constant. Second, the higher the terminal penalty $\varrho$ is, the relatively lesser weight is assigned to the generations up to time T, and therefore the retirement income level will be lower. Since $\bar{b}_L$ is unaffected by $\varrho$ , this implies a lower $\bar{b}_F$ and thus a higher indexation level. The same argument applies to the time preference rate since a lower $\zeta$ assigns relatively lower weights to the current generation. Due to the insignificant effects of $\zeta$ and $\varrho$ , we set $\zeta = 0$ and $\varrho = 1$ for the numerical analysis later on.

Figure 5. Partial indexation $\frac{\bar{b}_L\times L(t) }{\bar{b}_F + \bar{b}_L \times L(t) }$ for $t = 0$ .

We emphasize that the values of the TB level $\bar{b}(t)$ depend on the sponsors’ subjective choices of the risk-sharing thresholds $\xi_A$ and $\xi_{\text{VIX}}$ . A lower TB level is anticipated as a result of lower criteria because they imply that surplus payouts are triggered more frequently. There is no such thing as a free lunch, so it is unrealistic to anticipate having both a bigger surplus distribution and a higher TB level.

5.1. Salary uncertainty

Figure 6 presents the effects of the salary risk on the percentage of partial indexation and the volatility adjustment $\beta_{\text{VIX}}$ over different values of the mean-reverting speed $\kappa_l$ for the salary process and different values of the salary volatility $\sigma_l$ . It can be observed from the left panel that the mean-reverting speed $\kappa_l$ is crucial in determining the partial indexation level while the salary volatility is rather insignificant. We can also observe that the partial indexation stays within a practical range for different values of $\kappa_l$ and $\beta_{\text{VIX}}$ (i.e., approximately between 35% and 55%).

Figure 6. Impact of salary risk on partial indexation $\frac{\bar{b}_F}{\bar{b}_F+ \bar{b}_L }$ and $\beta_{\text{VIX}}$ .

The partial indexation level decreases as the mean-reverting speed $\kappa_l$ increases. This might be because the salary L is more likely to be close to the long-term mean $\overline{L}$ under a larger $\kappa_l$ and thus indicates that the benefit payment is less related to the wage index. The impact of the salary volatility $\sigma_l$ on the partial indexation is rather insignificant compared with $\kappa_l$ . In addition, the volatility adjustment $\beta_{\text{VIX}}$ is shown to be stable if the salary index is close to what is anticipated (when $\kappa_l$ is large or when $\sigma_l$ is small).

5.2. Jump risk

Figure 7 presents the effect of jump risks on $\beta_{\text{VIX}}(t)$ and $\bar{b}_F(t)$ via the jump intensity $\theta$ and the expected jump size $\mu_z$ . Overall, we can observe the monotone decreasing relationships between $\theta$ and both $\beta_{\text{VIX}}(t)$ and $\bar{b}_F(t)$ in the two panels, while the expected jump size $\mu_z$ has only insignificant impact. In the left panel, the higher jump risk premia allow the sponsor to set a higher overall TB level, and thus a higher $\bar{b}_F$ . Since $\bar{b}_L$ is unaffected, the level of partial indexation will be lower. In the right panel, the volatility adjustment is significantly affected by the amount of jump risk premia. Specifically, a higher $\theta$ represents both a higher frequency of downward shock and a higher jump risk premia, which in either case leads to a higher volatility adjustment. This observation aligns with our interpretation on the volatility adjustment in the previous section. Therefore, ignoring the relationship between the jump risk and the stochastic volatility would significantly underestimate the volatility adjustment in the optimal benefit structure.

Figure 7. Impact of jump risk on partial indexation $\frac{\bar{b}_F}{\bar{b}_F+ \bar{b}_L }$ and $\beta_{\text{VIX}}$ .

5.3. Dependency ratio

Most optimal control literature on pension plans assume a stationary population, which clearly does not match with the current experience of the aging society. We examine the impact of the aging population on both the performance adjustment and the volatility adjustment in this section.

Figure 8. United Nation Projection in old-age dependency ratio.

Based on the projection of United Nations (2022), the population growth for developed regions is virtually zero. Therefore, to study the impact of the changing population structure, we set $\mathcal{A}(t) = \mathcal{A}$ and solely focus on the improvement in the old age mortality rate. Figure 8 plots the projections of the old-age dependency ratioFootnote ¹⁸ for the selected developed countries. Although the projection values are significantly different across different countries, they generally are of an S-shape, which means most of the countries are experiencing a rapidly growing aging population and the growth will eventually slow down. It is therefore natural to use a Sigmoid functionFootnote ¹⁹ to approximate the dynamic of the dependency ratio. Define the old-age dependency ratio as $\mathcal{D}(t) = \frac{\mathcal{R}(t)}{\mathcal{A}(t)}$ with ODE

\begin{align*} {\mathop{}\!\mathrm{d}} \mathcal{D}(t) = \epsilon \times \left( \frac{\bar{\mathcal{D}} - \mathcal{D}(t) }{\bar{\mathcal{D}}} \right) \times \mathcal{D}(t) {\mathop{}\!\mathrm{d}} t,\end{align*}

where $\epsilon$ is the base growth rate and $\bar{\mathcal{D} }$ is the maximum dependency ratio can be achieved.

Figure 9 displays the effects of population growth rate on both the performance adjustment $\beta_A$ and the volatility adjustment $\beta_{\text{VIX}}$ , where $\bar{D} = 73\%$ is based on the UN’s projection on Canadian data. From the left panel, we can observe a higher growth rate (corresponding to a larger retiree population throughout time) leads to a lower performance adjustment at initial periods of time (i.e., approximately 15–20 years) since the initial population size is relatively small. As time t increases, when the growth of the retiree population slows down, the relative position of the population at time t may become smaller for large $\epsilon$ and the population experiences a larger performance adjustment. Notice also that as long as $\epsilon>0$ , the difference is rather insignificant in the sense that $\beta_A$ is rather stable once the long-term projection of the dependency ratio can be accurately estimated.

Figure 9. United Nation Projection in old-age dependency ratio.

The right panel displays a significantly different trend in volatility adjustment under different population dynamics, compared to the performance adjustment. However, if we scale the volatility adjustment by the number of retirees $\left(\frac{\beta_{\text{VIX}}(t)}{\mathcal{R}(t) } \right)$ to examine the volatility adjustment for each individual, then the difference is much less significant. This is not surprising since the volatility index $\text{VIX}_t^2$ only appears in measuring the utility of each retirement benefit not in the terminal penalty. Thus, $\beta_A$ is rather stable in terms of the aggregate adjustment, while $\beta_{\text{VIX}}$ is more stable in terms of the individual adjustment.