2.1. GP regression

Given a dataset $\mathcal{D} = \{\textbf{x}_i, y_i\}_{i=1}^n$ , the GP assumption and observation likelihood (2.3) imply $[\textbf{f}, \textbf{y}]^\top \sim \mathcal{MVN}$ , where $\textbf{f} = [f(\textbf{x}_1), \ldots, f(\textbf{x}_n)]^\top$ and $\textbf{y} = [y(\textbf{x}_1), \ldots, y(\textbf{x}_n)]^\top$ , so that the posterior $\textbf{f} | \textbf{y} \sim \mathcal{MVN}$ as well. More generally, for $\textbf{x}, \textbf{x}^{\prime} \in \mathbb{R}^d$ , $[f(\textbf{x}), f(\textbf{x}^{\prime}), \textbf{y}]^\top \sim \mathcal{MVN}$ , so that the posterior finite dimensional distribution is fully known:

(2.5) \begin{equation} [f_*(\textbf{x}), f_*(\textbf{x}^{\prime})]^\top \,:\!=\, \left([f(\textbf{x}), f(\textbf{x}^{\prime})]^\top | \textbf{y}\right) \sim \mathcal{MVN} \!\left([m_*(\textbf{x}), m_*(\textbf{x}^{\prime})]^\top, \begin{bmatrix} k_*(\textbf{x}, \textbf{x}) & \quad k_*(\textbf{x}, \textbf{x}^{\prime}) \\[5pt] k_*(\textbf{x}^{\prime}, \textbf{x}) & \quad k_*(\textbf{x}^{\prime}, \textbf{x}^{\prime}) \end{bmatrix}\right), \end{equation}

where, for arbitrary $\textbf{x}, \textbf{x}^{\prime} \in \mathbb{R}^d$ , the posterior covariance kernel is defined as $k_*(\textbf{x}, \textbf{x}^{\prime}) \,:\!=\, \textrm{cov}(f(\textbf{x}), f(\textbf{x}^{\prime}) | \textbf{y})$ . The Kolmogorov extension theorem ensures that $\{f_*(\textbf{x})\}_{\textbf{x} \in \mathbb{R}^d}$ defines a GP. Furthermore, (setting $m(\textbf{x}) \equiv 0$ temporarily) the posterior mean and variance are explicitly given by

(2.6) \begin{align} m_*(\textbf{x}) &= \textbf{K}(\textbf{x}, \textbf{X})\left[\textbf{K}(\textbf{X}, \textbf{X}) + \Delta(\textbf{X}, \textbf{X})\right]^{-1}\textbf{y}, \end{align}

(2.7) \begin{align} \textbf{K}_*(\textbf{x}, \textbf{x}^{\prime}) &= \textbf{K}\!\left([\textbf{x}, \textbf{x}^{\prime}]^\top, \textbf{X}\right) \left[\textbf{K}(\textbf{X}, \textbf{X}) + {\Delta}(\textbf{X}, \textbf{X})\right]^{-1} \textbf{K}\!\left(\textbf{X}, [\textbf{x}, \textbf{x}^{\prime}]^\top\right),\end{align}

where $\textbf{X}$ denotes the $n \times d$ matrix with rows $\textbf{x}_i, i=1, \ldots, n$ , and for $\textbf{U}, \textbf{V}$ being $\ell \times d$ and $m \times d,$ respectively, $\textbf{K}(\textbf{U}, \textbf{V})=\left[k(\textbf{u}_i, \textbf{v}_j)\right]_{1 \leq i \leq \ell, 1 \leq j \leq m}$ denotes the $\ell \times m$ matrix of pairwise covariances. ${\Delta}(\textbf{U}, \textbf{V})$ has entries $\sigma^2(\textbf{u}_i) \delta_{i = j}\delta_{\textbf{u}_i = \textbf{v}_j}$ ; in our case ${\Delta}(\textbf{X}, \textbf{X})$ is a $n \times n$ diagonal matrix with entries $\sigma^2(\textbf{x}_i)$ . In the case of constant noise variance $\sigma^2 \,:\!=\, \sigma^2(\textbf{x}_i)$ , ${\Delta}(\textbf{X}, \textbf{X}) = \sigma^2 I_n$ , where $I_n$ is the $n \times n$ identity matrix.

2.2. GP kernels

GP regression encodes the idea that similar inputs (according to the kernel) yield similar outputs. This can be seen through the posterior mean being a weighted average of observed data, since $m_*(\textbf{x}) = \textbf{w}^\top \textbf{y}$ holds for $\textbf{w}^\top = \textbf{K}(\textbf{x}, \textbf{X})\left[\textbf{K}(\textbf{X}, \textbf{X}) + \Delta(\textbf{X},\textbf{X})\right]^{-1}$ . Various types of kernels exist to encode similarity according to domain knowledge. Akin to covariance matrices, the only requirement for a function $k({\cdot}, \cdot)$ to be a covariance kernel is that it is symmetric and positive-definite, that is for all $n = 1, 2, \ldots,$ and $x_1, \ldots, x_n \in \mathbb{R}^d$ , we must have the Gram matrix $\textbf{K}(\textbf{X}, \textbf{X})$ be positive semi-definite.

A stationary kernel $k\big(\textbf{x}_i, \textbf{x}_j\big)$ is one that can be written as a function of $\textbf{x}_i - \textbf{x}_j$ , that is, $k\big(\textbf{x}_i, \textbf{x}_j\big) = k_S(\textbf{x}_i - \textbf{x}_j) $ and is thus invariant to translations in the input space. A kernel is further called isotropic if it is only a function of $\textbf{r} = \Vert\textbf{x}_i - \textbf{x}_j\Vert$ , where $\Vert \cdot \Vert$ is the $\ell^2$ Euclidean distance, so that we can write $k_I(\textbf{r}) = k\big(\textbf{x}_i, \textbf{x}_j\big)$ . Yaglom (Reference Yaglom1957) uses Bochner’s theorem to derive a similar Fourier transform specific to isotropic kernels, providing a way to derive kernels from spectral densities. Stationary kernels are usually assumed to be normalized, since $k_S(\textbf{x}-\textbf{x}^{\prime})/k_S(\textbf{0})=1$ whenever $\textbf{x} = \textbf{x}^{\prime}$ ; this allows for stationary kernels to be nicely interpreted as correlation functions. Let $\sigma_f^2$ denote the process variance $\textrm{var}(f(\textbf{x})) = \sigma_f^2$ , we think of the GP f as $f(\textbf{x}) = \sigma_f g(\textbf{x})$ when g is a GP with normalized stationary kernel k.

Lastly, a separable kernel over $\mathbb{R}^d$ is one that can be written as a product: $k(\textbf{x}, \textbf{x}^{\prime}) = \prod_{j=1}^d k_j\!\left(x^{(j)}, x^{^{\prime}(j)}\right)$ , where $x^{(j)}$ is the jth coordinate of $\textbf{x}$ . Thus, the global kernel is separated as a product over its dimensions, each having its own kernel.

The algebraic properties of positive-definite functions make it straightforward to compose new kernels from existing ones (Genton, Reference Genton2001; Schölkopf et al., Reference Schölkopf, Smola and Bach2002; Shawe-Taylor and Cristianini, Reference Shawe-Taylor and Cristianini2004; Berlinet and Thomas-Agnan, Reference Berlinet and Thomas-Agnan2011). The main tool is that kernels are preserved under addition and multiplication, that is, can be combined by sums and products. Hence, if $k_1$ and $k_2$ are two kernels and $c_1, c_2$ are two positive real numbers, then so is $k(\textbf{x}, \textbf{x}^{\prime}) = c_1 k_1(\textbf{x}, \textbf{x}^{\prime}) + c_2 k_2(\textbf{x}, \textbf{x}^{\prime})$ . This is consistent with the properties of GPs: if $f_1\sim \mathcal{G}\mathcal{P}(0, k_1)$ and $f_2 \sim \mathcal{G}\mathcal{P}(0, k_2)$ are two independent GPs, then for $c_1, c_2 >0$ we have $c_1 f_1+c_2 f_2 \sim \mathcal{G}\mathcal{P}\!\left(0, c_1 k_1 + c_2 k_2\right)$ . This offers a connection to the framework of generalized additive models. Using this, a GP whose kernel function is of the form $\sigma_{f_1}^2 k_1 + \sigma_{f_2}^2 k_2$ where $k_1$ and $k_2$ are normalized stationary kernels can be interpreted as $f(\textbf{x}) = \sigma_{f_1} f_1(\textbf{x}) + \sigma_{f_2} f_2(\textbf{x})$ where $f_1$ and $f_2$ are independent GPs with respective kernels $k_1$ and $k_2$ . Consequently, $\textrm{var}(f(\textbf{x})) = \sigma_{f_1}^2 + \sigma_{f_2}^2$ . These properties extend inductively to the case of finitely many terms.

Although there is no analogous result for a product of kernels (a product of GPs is no longer a GP, since the multivariate Gaussian distribution is not preserved), a common interpretation of multiplying kernels occurs when one kernel is stationary and monotonically decaying as $|\textbf{x}-\textbf{x}^{\prime}| \rightarrow \infty$ . Indeed, if $k_1$ is such a kernel, then the product $k_1(\textbf{x},\textbf{x}^{\prime})k_2(\textbf{x},\textbf{x}^{\prime})$ offers $k_2$ ’s effect with a decay according to $k_1(\textbf{x},\textbf{x}^{\prime})$ as $|\textbf{x}-\textbf{x}^{\prime}|$ increases. This synergizes with separability, for example, $k(\textbf{x}, \textbf{x}^{\prime}) = k_1\!\left(|x^{(1)} - x^{^{\prime}(1)}|\right)k_2\!\left(x^{(2)}, x^{^{\prime}(2)}\right)$ offers some similarity across the second coordinate that could decay as $|x^{(1)} - x^{^{\prime}(1)}|$ increases. Additive and multiplicative properties are often used in conjunction with the constant kernel $k(\textbf{x}, \textbf{x}^{\prime}) = c, c >0$ , resulting in a scaling effect (multiplication) or dampening effect (addition).

2.3. Kernel families

Table 1 lists the kernel families we consider. For simplicity, we assume one dimensional base kernels with $x \in \mathbb{R}$ where the full structure for $\textbf{x} \in \mathbb{R}^d$ is expressed as a separable kernel. Among the nine kernel families, we have kernels that give smooth ( $C^2$ and higher) fits, kernels with rough (non-differentiable) sample paths, and several nonstationary kernels. Here, smoothness refers to the property of the sample paths $\textbf{x} \mapsto f(\textbf{x})$ being, say, k-times differentiable, that is, $f({\cdot}) \in C^k$ . The posterior mean $\textbf{x} \mapsto m_*(\textbf{x})$ inherits similar (but generally less strict) differentiability properties, see Kanagawa et al. (Reference Kanagawa, Hennig, Sejdinovic and Sriperumbudur2018) for details. The $\mathcal{K}_r$ column indicates whether a kernel is included in the restricted set of kernels (used in later sections), in contrast to the full set of kernels $\mathcal{K}_f$ . This subset $\mathcal{K}_r$ comprises a more compact collection of the most commonly used kernels in the literature. One reason for $\mathcal{K}_r$ is to minimize overlap in terms of kernel properties; as we report below, some kernel families in $\mathcal{K}_f$ apparently yield very similar fits and act as “substitutes” for each other.

Table 1. List of kernel families used in compositional search. $C^p$ indicates that the GP sample paths $x \mapsto f(x)$ have p continuous derivatives; $C^0$ is continuous but not differentiable. Column $\mathcal{K}_r$ denotes whether the kernel family is in the restricted search set. The linear kernel is used for its year component only.

All of the kernels listed have hyperparameters, which help to understand their relationship with the data. The quantity $\ell_{\textrm{len}}$ appearing in many stationary kernels is referred to as the characteristic lengthscale, which is a distance-scaling factor. With the radial basis function (RBF) kernel for example, this loosely describes how far x ^′ needs to move from x in the input space for the function values to become uncorrelated (Williams and Rasmussen, Reference Williams and Rasmussen2006). Thus, the lengthscale of a calibrated GP can be interpreted as the strength of the correlation decay in the training dataset. Out of the stationary kernels, a popular class is the Matérn class. In continuous input space, the value $\nu$ in the Matérn- $\nu$ corresponds to smoothness: a GP with a Matérn- $\nu$ kernel is $\lceil \nu \rceil -1$ times differentiable in the mean-square sense (Williams and Rasmussen, Reference Williams and Rasmussen2006). The RBF kernel is the limiting case as $\nu \rightarrow \infty$ , resulting in an infinitely differentiable process. The $\nu=1/2$ case recovers the well-known Ornstein–Uhlenbeck process, which is mean reverting and non-differentiable. Also non-differentiable but on the nonstationary side, the minimum kernel corresponds to a Brownian motion process when $x \in \mathbb{R}_+$ , where $t_0^2 = \textrm{var}(f(0))$ is the initial variance. For discrete x, Min kernel yields the random walk process, and $\textrm{M}12$ yields the AR(1) process.

Less commonly studied are the (continuous-time) AR2 (see Parzen, Reference Parzen1961), Cauchy, and Mehler kernels. The continuous-time AR2 kernel acts identically to a discrete-time autoregressive (AR) process of order 2 with complex characteristic polynomial roots when x is restricted to an integer. The heavy-tailed Cauchy probability density function motivates the Cauchy kernel, with the goal of modeling long-range dependence and is a special case of the rational quadratic kernel (see Appendix A). Lastly, the Mehler kernel has a form similar to that of a joint-normal density and acts as a RBF kernel with a nonstationary modification (this can be seen from a “complete the squares” argument). Although Mehler is nonstationary, it remarkably yields a stationary correlation function $\textrm{corr}(f(\textbf{x}), f(\textbf{x}^{\prime})) = k(\textbf{x}, \textbf{x}^{\prime}) / \sqrt{k(\textbf{x}, \textbf{x}) k(\textbf{x}^{\prime}, \textbf{x}^{\prime})}$ . See Appendix A for a more thorough discussion of the aforementioned kernels and their properties.

2.4. Varying prior means

For brevity of exposition, the preceding analysis considered zero prior mean $m(\textbf{x}) = 0$ . When the prior mean function $m(\textbf{x})$ is known and deterministic, the posterior covariance function $k_*(\textbf{x}, \textbf{x}^{\prime})$ remains unchanged compared to the zero-mean case, and the posterior mean $m_*(\textbf{x})$ is adjusted to:

(2.8) \begin{equation}m_*(\textbf{x}) = m(\textbf{x}) + \textbf{K}(\textbf{x}, \textbf{X})\left[\textbf{K}(\textbf{X}, \textbf{X}) + \Delta(\textbf{X}, \textbf{X})\right]^{-1}(\textbf{y} - m(\textbf{X})),\end{equation}

where $m(\textbf{X}) = [m(\textbf{x}_1), \ldots, m(\textbf{x}_n)]^\top$ . More common is the case of a parametric mean function $m(\textbf{x}) = \sum_{j=1}^d \beta_j h_j(\textbf{x})$ , where the basis functions $h_j(\textbf{x})$ are fixed and known (e.g., $h_j(\textbf{x})$ is a jth degree polynomial), and the coefficients $\boldsymbol{\beta}$ are estimated simultaneously with the covariance hyperparameters through maximizing the marginal likelihood. The generalized least squares estimator for $ \boldsymbol{\beta} $ is

(2.9) \begin{equation}\hat{\boldsymbol{\beta}} = \left( \textbf{H}^\top \textbf{K}_{\textbf{y}}^{-1} \textbf{H} \right)^{-1} \textbf{H}^\top \textbf{K}_{\textbf{y}}^{-1} \textbf{y}, \qquad \textrm{with} \quad \textbf{K}_{\textbf{y}} = \textbf{K}(\textbf{X}, \textbf{X}) + \Delta(\textbf{X}, \textbf{X}),\end{equation}

where $\textbf{H}$ is the $n \times d$ matrix with i, j entry $h_j(\textbf{x}_i)$ . The resulting posterior mean $ m_*(\textbf{x}) $ is based on the plug-in estimated trend $\textbf{H}\hat{\boldsymbol{\beta}}$ and the posterior variance has an additional term to account for the parameter uncertainty in $\boldsymbol{\beta}$ :

\begin{align*}m_*(\textbf{x};\, \hat{\beta}) &= \textbf{h}(\textbf{x})^\top \hat{\boldsymbol{\beta}} + \textbf{K}(\textbf{x}, \textbf{X})\textbf{K}_{\textbf{y}}^{-1}\left(\textbf{y} - \textbf{H}\hat{\boldsymbol{\beta}}\right), \qquad \textbf{h}(\textbf{x}) = [h_1(\textbf{x}), \ldots, h_d(\textbf{x})]^\top;\\k_*(\textbf{x}, \textbf{x};\, \hat{\beta}) &= k_*(\textbf{x}, \textbf{x}) +\left(\textbf{h}(\textbf{x}) - \textbf{K}(\textbf{x}, \textbf{X})^\top \textbf{K}_{\textbf{y}}^{-1} \textbf{H}\right)^\top\left(\textbf{H}^\top \textbf{K}_{\textbf{y}} \textbf{H}\right)^{-1}\left(\textbf{h}(\textbf{x}) - \textbf{K}(\textbf{x}, \textbf{X})^\top \textbf{K}_{\textbf{y}}^{-1} \textbf{H}\right).\end{align*}

For more details, see Roustant et al. (Reference Roustant, Ginsbourger and Deville2012). It is important to recall that $\textbf{x}$ is an input to the function, while $\textbf{X}$ is fixed (from the observed data). Intuitively, this reflects a transition toward prior reliance: in the case of a decaying kernel (e.g., Matérn family), when the predictive location $\textbf{x}$ distances itself from the rows of $\textbf{X}$ , that is, $\Vert\textbf{x} - \textbf{x}_i\Vert \rightarrow 0$ for all training locations $\textbf{x}_i$ , $i = 1, \ldots, n$ , the vector $\textbf{K}(\textbf{x}, \textbf{X})=[k(\textbf{x}, \textbf{x}_i)]_{i=1}^n$ goes to $\textbf{0}$ , causing a reversion to the prior mean and covariance.

2.5. Connections to mortality modeling

For mortality modeling, our core input space is composed of Age and Year coordinates: $x_{ag}, x_{yr} \in \mathbb{R}_+^2$ . As the GP can model non-linear relationships, we include Birth Cohort $x_c \,:\!=\, x_{ag}-x_{yr}$ as a third coordinate of $\textbf{x}$ , so that $\textbf{x} = (x_{ag}, x_{yr}, x_{c})$ . For a given $\textbf{x}$ , denote $D_{\textbf{x}}$ and $E_{\textbf{x}}$ as the respective observed deaths and exposures (i.e., individuals alive at the beginning of the period) over the corresponding $(x_{ag}, x_{yr})$ pair. Denote $y(\textbf{x}) = \log\!(D_{\textbf{x}} / E_{\textbf{x}})$ as the log mortality rate. The full dataset is denoted $\mathcal{D} = \{\textbf{x}_i, y_i, D_i, E_i\}_{i=1}^n$ . Each mortality observation $y_i \,:\!=\, y_i(\textbf{x}_i)$ is the regression quantity modeled in Equation (2.3), so that $\mathbb{E}[y(\textbf{x}) | f(\textbf{x})] = f(\textbf{x})$ . The interpretation is that the true mortality rate $f(\textbf{x})$ is observed in accordance with mean-zero (i.e., unbiased) uncorrelated noise yielding the measured mortality experience y.

Relating to the log-normal distribution, algebra shows that

\begin{equation*}\mathbb{E}[D_{\textbf{x}} | f(\textbf{x})] = E_{\textbf{x}} \exp\!\left(f(\textbf{x}) + \frac{\sigma^2(\textbf{x})}{2}\right), \textrm{ and } \quad \mathbb{E}[D_{\textbf{x}} | f(\textbf{x}), \epsilon(\textbf{x})] = E_{\textbf{x}}\exp\!\left(f(\textbf{x}) + \epsilon(\textbf{x})\right),\end{equation*}

akin to an overdispersed Poisson model in existing mortality modeling literature (see, e.g., Azman and Pathmanathan, Reference Azman and Pathmanathan2022). Note that the seminal work of Brouhns et al. (Reference Brouhns, Denuit and Vermunt2002) for modeling log-mortality rates suggests that homoskedastic noise is unrealistic, since the absolute number of deaths at older ages is much smaller compared to younger ages. As a result, we work with heteroskedastic noise

(2.10) \begin{align}\textrm{var}(\epsilon | D_{\textbf{x}}) = \sigma^2(\textbf{x}) \,:\!=\, \frac{\sigma^2}{D_{\textbf{x}}}, \quad \textrm{where } \sigma^2 \in \mathbb{R}^+.\end{align}

Thus, we make observation variance inversely proportional to observed death counts, with the constant $\sigma^2$ to be learned as part of the fitting procedure. From a modeling perspective, this works since $D_{\textbf{x}}$ is known whenever $E_{\textbf{x}}$ and $y(\textbf{x})$ are. Indeed, Equation (2.5) shows no requirement to know $D_{\textbf{x}}$ for out-of-sample forecasting. In the case where full distributional forecasts are desired, one could instead model a noise surface $\sigma^2({\cdot})$ simultaneously with $f({\cdot})$ , see, for example, Cole et al. (Reference Cole, Gramacy and Ludkovski2022).

A discussion of GP covariances connects naturally to APC models. A stationary covariance means that the dependence between different age groups or calendar years is only a function of the respective ages/year distances and is not subject to additional structural shifts. Additive and separable covariance structures play an important role in the existing mortality modeling literature, specifically in the APC framework. For example, Lee–Carter and CBD began with Age–Period modeling, which subsequently evolved to add cohort or additional Period effects. Rather than postulating the precise APC terms, the latest (Dowd et al., Reference Dowd, Cairns and Blake2020) CBDX framework adds up to 3 period effects as needed. Similarly, Hunt and Blake (Reference Hunt and Blake2014) develop a general recipe for constructing mortality models, where core demographic features are represented with a particular parametric form, and combined into a global structure. This can be reproduced with GPs where kernels encode expert judgment to such demographic features. This type of encoding into covariances is already being done in the literature, possibly unknowingly. Take, for example, the basic CBD model whose stochastic part has the form $f(\textbf{x}) = \kappa_{x_{yr}}^{(1)} + \left(x_{ag}-\overline{\textbf{x}_{ag}}\right) \kappa_{x_{yr}}^{(2)}$ , where $\overline{\textbf{x}_{ag}}$ is the average age in $\mathcal{D}$ . Under the common assumption that $\left(\kappa_{x_{yr}}^{(1)}, \kappa_{x_{yr}}^{(2)}\right)$ is a multivariate random walk with drift, a routine calculation shows that

\begin{equation*}\mathbb{E}[f(\textbf{x})] = \kappa_0^{(1)}+\mu^{(1)}x_{yr} + \left(x_{ag}-\overline{\textbf{x}_{ag}}\right)\left(\kappa_0^{(2)}+\mu^{(2)}x_{yr}\right)\end{equation*}

and $ k_{CBD}(\textbf{x}, \textbf{x}^{\prime}) = \textrm{cov}(f(\textbf{x}), f(\textbf{x}^{\prime})) =$

(2.11) \begin{equation} = \left[ \sigma_1^2 + \rho \sigma_1 \sigma_2\!\left(x_{ag} + x^{\prime}_{ag} - 2\overline{\textbf{x}_{ag}}\right) + \left(x_{ag} -\overline{\textbf{x}_{ag}}\right)\left(x^{\prime}_{ag} - \overline{\textbf{x}_{ag}}\right) \sigma_2^2 \right]\left(x_{yr} \wedge x^{\prime}_{yr}\right).\end{equation}

We can re-interpret the above as the kernel decomposition $k_{CBD}(x, x^{\prime}) = k_1\!\left(x_{ag}, x^{\prime}_{ag}\right)k_2\!\left(x_{yr}, x^{\prime}_{yr}\right)$ , where $k_2\!\left(x_{ag}, x^{\prime}_{ag}\right)$ is a kernel depending only on Age (that could be decomposed into additive components), and $k_2$ is the minimum kernel depending only on Year. Furthermore, if the multivariate random walk is assumed to be Gaussian, then $f(\textbf{x})$ actually forms a (discrete-time) GP, and hence the methods detailed in Section 2.1 apply verbatim.

Period and Cohort effects are commonly modeled using time series models (Villegas et al., Reference Villegas, Kaishev and Millossovich2015). In particular, Gaussian ARIMA models are popular in existing APC literature, see, for instance, Cairns et al. (Reference Cairns, Blake, Dowd, Coughlan, Epstein and Khalaf-Allah2011) who single out the usefulness of AR(1), ARIMA(1,1,0), and ARIMA(0,2,1) for Cohort effect, and ARIMA(1,1,0) or ARIMA(2,1,0) for Period effect. This provides another link to (discrete) GP covariance analogues: a Matérn-1/2 covariance corresponds to an AR(1) process, a (continuous-time) AR2 covariance to an AR(2) process with complex unit roots, and a Minimum covariance to a discrete-time random walk.

3.1. BIC and Bayes factors for GPs

To evaluate the appropriateness of a kernel within a given set $k \in \mathcal{K}$ , an attractive criterion is the posterior likelihood of the kernel given the data $p(k | \textbf{y}) = p(\textbf{y} | k) p(\textbf{y}) / p(k)$ , where, under a uniform prior assumption $p(k) = 1/|\mathcal{K}|$ , we see that $p(k | \textbf{y}) \propto p(\textbf{y} | k)$ . However, the integral over hyperparameters $p(\textbf{y} | k) = \int_\theta p(\textbf{y}, \theta| k)d\theta$ is generally intractable, so we use the BIC as an approximation, where $\textrm{BIC}(k) \approx \log p(\textbf{y} | k)$ is defined as:

(3.1) \begin{equation} \textrm{BIC}(k) = -l_{k}\left(\hat{\theta};\, \textbf{y}\right) + \frac{|\hat{\theta}| \log\!(n)}{2},\end{equation}

where $l_k(\theta| \textbf{y}) = \log p(\textbf{y} | k, \theta)$ is the log marginal likelihood of y evaluated at $\theta$ under a given kernel k, $\hat{\theta}$ is the maximizer (maximum marginal likelihood estimate) of $l_k(\theta| \textbf{y})$ , and $|\hat{\theta}|$ is the total number of estimated hyperparameters in $\hat{\theta} = \left[\hat{\beta}_0, \hat{\beta}_{ag}, \hat{\theta}_k, \hat{\sigma}^2\right]^\top$ , where $\theta_k$ is a vector of all kernel specific hyperparameters. Note that $p(\textbf{y}|k, \theta)$ is a multivariate density, with mean and covariances governed by Equation (2.4). The BIC metric has seen use in similar applications of GP compositional kernel search, see, for example, Duvenaud et al. (Reference Duvenaud, Lloyd, Grosse, Tenenbaum and Zoubin2013), Duvenaud (Reference Duvenaud2014), and is commonly used in mortality modeling (Cairns et al., Reference Cairns, Blake, Dowd, Coughlan, Epstein, Ong and Balevich2009). We employ BIC (lower BIC being better) for our GA fitness metric below.

One can further assess the relative likelihood of $k_1, k_2 \in \mathcal{K}$ by again assuming a uniform prior over $\mathcal{K}$ , and computing the Bayes factor (BF)

(3.2) \begin{equation} \textrm{BF}(k_1, k_2) = \frac{p(k_1 | \textbf{y})}{p(k_2 | \textbf{y})} \approx \exp\Big(\textrm{BIC}(k_2) - \textrm{BIC}(k_1)\Big).\end{equation}

Obtaining an approximation for $p(k | \textbf{y})$ further allows the possibility of Bayesian model averaging by conditioning over $k \in \mathcal{K}$ :

(3.3) \begin{equation} p(\textbf{f} | \textbf{y}) = \sum_{k \in \mathcal{K}} p(\textbf{f} | \textbf{y}, k) p(k | \textbf{y}).\end{equation}

Since $p(\textbf{f} | \textbf{y}, k)$ is available in the closed form, this provides the full distribution of future forecasts, if desired.

Figure 1. Representative compositional kernels and GA operations. Bolded red ellipses indicate the node of $\kappa$ (or $\xi$ ) that was chosen for mutation or crossover.

Gelman et al. (Reference Gelman, Carlin, Stern and Rubin1995) states that BFs work well in the case of a discrete model selection. The seminal work of Jeffreys (Reference Jeffreys1961) gives a table of evidence categories to determine a conclusion from BFs, see Table 15 in Appendix E, which is still frequently used today (Lee and Wagenmakers, Reference Lee and Wagenmakers2014; Dittrich et al., Reference Dittrich, Leenders and Mulder2019).

Note that in accordance with the penalty term $|\hat{\theta}| \log\!(n)/2$ in Equation (3.1), the difference in penalties in $\textrm{BIC}(k_2) - \textrm{BIC}(k_1)$ is simply the number of additional kernel hyperparameters, scaled by $\log\!(n)/2$ , since $\beta_0, \beta_{ag}, \sigma^2$ are always estimated regardless of kernel choice. Thus in the application of GP kernel selection, the BF properly penalizes kernel complexity. In Table 1, most kernels have one $\theta_k$ (the lengthscale), but some, like $\textrm{AR}2$ , have two hyperparameters and so incur a larger penalty.

3.2. GA kernel representation

In order to operate in the space of kernels, we shall represent kernels via a tree-like structure, cf. Figure 1. Internal nodes correspond to operators (add or mul) that combine two different kernels together, while leafs are the univariate kernels used as building blocks. We indicate the coordinate operated on by the kernel through the respective subscripts a, y, c, such as $\textrm{M}52_{a}$ . In turn, such trees are transcribed into bracketed expressions, such as $\kappa = \texttt{add(Exp_c, mul(RBF_a, add(Mat_y, RBF_c)))}$ corresponding to the Age–Period–Cohort kernel $(k_{M52}(x_{yr})+k_{RBF}(x_c) )\cdot k_{RBF}(x_{ag}) + k_{Exp}(x_c)$ . The length of $\kappa$ , denoted $|\kappa|$ , is its number of nodes. The above kernel tree has length $|\kappa| = 7$ , namely four base kernels combined with three operators. Observe that the tree structure is not unique, that is, some complex kernels can be permuted and expressed through different trees. In what follows, we will ignore this non-uniqueness.

A certain expertise is needed to convert from an above representation to the dependence structure it implies. One way to visualize is to take advantage of the stationarity and plot the heatmap of the matrix $k_S(\textbf{x}-\textbf{x}^{\prime})$ as a function of its Age, Year coordinates.

3.3. GA operations and algorithm parameters

The overall GA is summarized by the set of possible mutations and a collection of algorithmic parameters. These are important for many reasons: (i) sufficient exploration so that the algorithm does not get trapped in a particular kernel configuration; (ii) efficiency in terms of number of generation and generation size needed to find the best kernels; and (iii) bloat control, that is, ensuring that returned kernels are not overly complex and retain interpretability. Interpretable kernels would tend to have low length (below 10) and avoid repetitive patterns. Bloat control, that is, avoiding the appearance of overly complex/long kernels, is a concern with GAs. Luke and Panait (Reference Luke and Panait2006) and Poli et al. (Reference Poli, Langdon, McPhee and Koza2008) suggest several ways to combat it.

Our specific high-level GA parameters are listed in Table 2. We largely follow guidelines from Sipper et al. (Reference Sipper, Fu, Ahuja and Moore2018), which offers a thorough investigation of the parameter space of GA algorithms. Ancestors are chosen via a tournament setup, where $T=7$ individuals are independently and uniformly sampled from the previous generation and a single tournament winner is the fittest (lowest BIC) individual. A smaller T provides diversity in future generations, whereas a larger T reduces chances of leaving behind fit individuals. To combat bloat, we follow the double tournament procedure described in Luke and Panait (Reference Luke and Panait2006): all instances of a single tournament are replaced by two tournaments run one after the other, with the lower length ancestor chosen with probability ${p_{DT} \in [0.5, 1]}$ ; larger values prefer parsimony over fitness. We use the less restrictive ${p_{DT}=0.6}$ instead of the suggested ${p_{DT}=0.7}$ , partially since BIC has a built-in penalization for unwieldy models thereby mitigating bloat. As a further proponent of parsimony, we use hoist mutation as suggested in Poli et al. (Reference Poli, Langdon, McPhee and Koza2008) with $p_h=0.1$ . Using the above parameterization, we rarely observe kernels of length over 15 in our experiments.

Table 2. High level Genetic algorithm parameters and description. Note that $G \cdot n_g = 4000$ is the total number of trained GP models in a single run of the GA.

Table 3 fully details the crossover and mutation operations and associated algorithmic parameters, with Figure 1 providing a visual illustration. Most of these are standard in the literature. Our domain knowledge suggests an additional point mutation operator which we call respectful point mutation. This operator maintains the coordinate of the kernel being mutated, so that a kernel operating on age is replaced by one in age, and so forth. This respects a discovered APC structure and fine tunes the chosen coordinate.

Table 3. Operator specific GA parameters with description. Note that $p_c + p_s + p_h + p_p + p_r + p_0 = 1$ .

The arity (number of arguments) of a function needs to be preserved in crossover and mutation operations. In our setup, the only non-trivial arity functions are add and mul (both with arity 2), so if these nodes are chosen for a mutation, the point mutation (and respectful point mutation) automatically replaces them with the other operation: add with mul and vice-versa.

Algorithm 1: Genetic Algorithm for compositional kernel selection

In some GA applications, authors propose to have several hundred (or even thousand) generations. Because times to fit a GP model is non-trivial, running so many generations is computationally prohibitive. Below we use a fixed number of $G=20$ generations. As shown in Figure 5, kernel exploration seems to stabilize after a dozen or so generations, so there appears to be limited gain in increasing G. In contrast with a typical GA application, we are more interested in interpreting all prototypical well-performing kernels, rather than in optimizing an objective function to an absolute minimum; that is to say, our analysis is sufficient as long as a representative ballpark of optimal models has been discovered.

Denote by $\kappa_i^g$ the $i = 1, \ldots, n_g$ th individual in generation $g=1, \ldots, G-1$ , with a corresponding BIC of $b_i^g$ . Algorithm 1 outlines the full GA. BIC computes the BIC as in Equation (3.1), which implicitly fits a GP and optimizes hyperparameters. SampleUniform determines ancestors according to a tournament of size T. A crossover or mutation is determined according to the probabilities provided in $\mathcal{G}_2$ (according to Table 3), where Type denotes the type of mutation chosen. In the event of a crossover, another ancestor is determined through the same process as the first. The entire algorithm is ran for G generations, where $g=0$ initializes, and the remaining $G-1$ steps involve choosing ancestors and offspring. For $g=0$ , InitializeKernel() constructs a randomly initialized kernel, where the length is uniformly chosen from $\{3, 5, 7, 9\}$ (respectively 2, 3, 4, 5 base kernels), where the base kernels are uniformly sampled from $\mathcal{K}$ , and add/multiply operations are equally likely. Experiments where we initialized with a more diverse set (kernel length from 1 to 15) slowed GA convergence with a negligible effect on end results. Our initialization was chosen with the goal of providing as unbiased of a sampling procedure as possible. Alternative initializations could be used, for example, imposing a diversification constraint on the APC terms in a given kernel, or infusing the initial generation with known-to-be-adequate kernels. It is also worth mentioning that although we fixed $n_g=200$ for all $g = 0, \ldots, G-1$ , this value could vary over g so that the total number of results is $\sum_{g=0}^{G-1} n_g$ .

Since top-performing kernels are preferred as tournament winners, they become ancestors for future generators and are likely to re-appear as duplicates (either due to a Copy operation or a couple of mutations canceling each other). As a result, we observe many duplicates when aggregating all $\kappa_i^g$ ’s across generations.

3.4. GP hyperparameter optimization

As mentioned, the hyperparameters of a given kernel $k({\cdot}, \cdot; \theta)$ are estimated by maximizing $l_k(\theta| \textbf{y})$ , the marginal log likelihood of the observed data. The optimization landscape of kernel hyperparameters in a GP is non-convex with many local minima, so we take care in this optimization. Since $\textbf{y}$ is on the log scale, we leave these values untransformed. For given $\textbf{x} = [x_{ag}, x_{yr}, x_{c}]^\top \in \mathbb{R}^3$ , we perform dimension-wise scaling to the unit interval, for example, if $x_{ag} = [x_{1,ag}, \ldots, x_{n,ag}]^\top$ then $x_{i,ag} \mapsto \frac{x_{i,ag} - \min(x_{ag})}{\max(x_{ag})-\min(x_{ag})}$ and similarly for $x_{yr}$ and $x_c$ . Lengthscales for stationary kernels can be interpreted on the original scale through the inverse transformation $\ell_{\textrm{len}} = \left(\max(x) - \min(x)\right) \tilde{\ell}_{\textrm{len}}$ , where $\tilde{\ell}_{\textrm{len}}$ is on the transformed scale, with similar transformations for the mean function parameters. For interpretability purposes, we utilize this to report values on the original scale in Sections 4 and 5 whenever possible. Nonstationary kernel (i.e., $\textrm{Min}, \textrm{Meh,}$ and Lin) hyperparameters are left transformed, as they generally cannot be interpreted on the original scale.

We use Python with the GPyTorch library (Gardner et al., Reference Gardner, Pleiss, Weinberger, Bindel and Wilson2018) to efficiently handle data and matrix operations. Since our results rely heavily on accurate likelihood values (through BIC), we turn off GPyTorch’s default matrix approximation methods. Maximizing $l_k(\theta| \textbf{y})$ is done using Adam (Kingma and Ba, Reference Kingma and Ba2014). This is an expensive procedure, as a naive evaluation of $l_k(\theta| \textbf{y})$ is $\mathcal{O}(n^3)$ from inverting $\textbf{K}(\textbf{X},\textbf{X})$ . Our GA runs use a convergence tolerance of $\varepsilon = 10^{-4}$ and maximum iterations of $\eta_{\textrm{max}} = 300$ . We found most simple kernels to converge quickly $(\eta \leq 100)$ even up to $\varepsilon = 10^{-6}$ , see Table B.1 in the Appendix. We keep $\eta_{\textrm{max}}$ relatively small since we need to fit $n_g \cdot G \gg 10^3$ models. Upon completion, the top few dozen kernels are refit with $\varepsilon = 10^{-6}$ and $\eta_{\textrm{max}} = 1000$ . Note that since Adam is a stochastic algorithm, the fitted kernel hyperparameters vary (slightly in our empirical work) during this refitting. This reflects the idea that the GA is a preliminary “bird’s-eye search” to find plausible mortality structures, thereafter refining hyperparameter estimates. Except in a few cases, changes in final BIC values are minimal, though the ranking of top kernels can occasionally get adjusted.

4.1. Synthetic results

Answers to questions presented in the previous section are found in Table 5 for SYA, Table 7 for SYB, and Table 8 for SYC. In all tables, $K_0$ denotes the known kernel that generated the synthetic data. With the goal in mind to establish kernel recovery (ignoring hyperparameter estimation), $\hat{\theta}$ is estimated for $K_0$ from the generated data.

Table 5. Top five fittest non-duplicate kernels for the first synthetic case study SYA. Bolded is $K_0 = \textrm{RBF}_y \textrm{RBF}_a$ , the true kernel used in data generation. SYA-1 and -2 denote the realization trained on.

SYA

Table 5 shows the results of the top 5 fittest kernels for SYA-1 and SYA-2. The estimated BF (using BIC) is the column $\widehat{\textrm{BF}}(k, K_0) = \exp(\textrm{BIC}(K_0) - \textrm{BIC}(k))$ . For SYA-1, the true kernel $K_0 = \textrm{RBF}_y \textrm{RBF}_a$ is discovered and appears with lowest BIC. Next best (with a large BF of $0.826$ ) is $\textrm{M}52_a \textrm{RBF}_y$ , which has an identical structure aside from using the slightly less smooth Matérn- $5/2$ kernel for age instead of RBF. Note that the BF would need to be below $1/3 \approx 0.3333$ to even be worth mentioning a difference between these kernels according to Table 15 in the Appendix, suggesting an indifference of the two models. Our interpretation is that sampling variability makes $M52_a$ a similar alternative to $RBF_a$ . Note, however, that any lower Matérn order does not appear in either table. The remaining three rows have a similar kernel structure with superfluous Age or Period kernels added. The BFs are below 0.1, suggesting strong evidence against these kernels as being plausible alternatives for the generated dataset. This experiment confirms the GA’s ability to recover the overall structure, with an understandable difficulty in distinguishing between $M52_a$ and $\textrm{RBF}_a$ . SYA-2 in the right panel of Table 5 considers a different realized mortality surface under the same population distribution to assess stability across GA runs. It tells a similar story, though interestingly the lowest-BIC kernel is now $\textrm{M}52_a \textrm{RBF}_y$ , with $\widehat{BF}(k, K_0) = 1.1907$ . This means that it achieves a lower BIC than the true kernel, showcasing possibility of overfitting to data. At the same time, since the BF is so close to 1, this is still not statistically worth mentioning and hence can be fully chalked up to sampling variability. Otherwise, we again observe only two truly plausible alternatives, and very similar less-plausible (BF between 0.05 and 0.13) alternates.

To further assess smoothness detection, multiplicative Age–Period kernels $k=k_a k_y$ are fit to both SYA-1 and SYA-2 datasets, where $k_a$ and $k_y$ range over $\textrm{M}12, \textrm{M}32, \textrm{M}52$ and RBF (in order of increasing smoothness), resulting in a total of 16 combinations per training surface. The resulting differences in BIC are provided in Table 6. Note that the ground truth kernel $K_0 = \textrm{RBF}_a \textrm{RBF}_y$ generates a mortality surface that is infinitely differentiable in both Age and Period. In both cases, decisive evidence is always against either surface having a $\textrm{M}12$ component (with $\textrm{BIC}(K_0) - \textrm{BIC}(k)$ ranging from $-96.99$ to $-30.24$ – recall that anything below –4.61 is decisive evidence for $K_0$ ). As found above, $\textrm{M}52_a$ is a reasonable surrogate for $\textrm{RBF}_a$ for both SYA-1 and SYA-2. SYA-1 shows only strong evidence against $\textrm{M}52_y$ (as opposed to decisive for SYA-2). In both cases, there is decisive evidence against $\textrm{M}32_y$ , with only strong evidence on the age component $\textrm{M}32_a$ when using $\textrm{RBF}_y$ (otherwise, it is decisive). This difference is likely explained by the additional flexibility of $\textrm{M}32_a$ to pick up some fluctuations which would normally be reserved for the parametric mean function in age.

Table 6. Logarithmic Bayes Factor $\log\widehat{BF}(k, K_0) = \textrm{BIC}(K_0) - \textrm{BIC}(k)$ for $K_0 = \textrm{RBF}_a \textrm{RBF}_y$ and $k = k_a k_y$ , where $k_a$ and $k_y$ are in the respective row and column labels and $\textrm{BIC}(K_0)$ is evaluated over SYA-1 (left panel) and SYA-2 (right panel).

Table 7. Top five fittest non-duplicate kernels for the second synthetic case study SYB. Bolded is the true kernel used in data generation $K_0 = 0.08 \cdot \textrm{RBF}_a({19.3}) \cdot M12_y({386.6}) + 0.02 \cdot M52_c({4.98})$ , with $\sigma^2 = 4 \times 10^{-4}, \beta_0 = {-9.942}, \beta_1 = {0.0875}$ .

Table 8. Fittest non-duplicate kernels for SYC in two separate runs, one over $\mathcal{K}_r$ and the other over $\mathcal{K}_f$ . The true kernel is $K_0 = \textrm{M}52_a \left(\textrm{Chy}_y \textrm{M}12_y\right) \textrm{M}12_c$ . Note that $\textrm{Chy}_y \notin \mathcal{K}_r$ . All $\widehat{\textrm{BF}}$ values are in the *** evidence category.

5.1. Robustness check

To validate the above results of the GA, we perform two checks: (i) re-run the algorithm from scratch, to validate stability across GA runs; (ii) run the GA on two modified training datasets: $\mathcal{D}_{\textrm{rob,1}}, \mathcal{D}_{\textrm{rob,2}}$ . For $\mathcal{D}_{\textrm{rob,1,}}$ we augment with two extra calendar years (beginning at 1988 instead of 1990), and four extra Age groups (48–86 instead of 50–84). For $\mathcal{D}_{\textrm{rob,2,}}$ we shift the dataset by 4 years in time, namely to 1986–2015, considering same Age range 50–84.

In all above cases, we expect results to be very similar to the “main” run discussed above. While the GA undertakes random permutations and has a random initialization, we expect that with 200 kernels per generation and 20 generations, the GA explores sufficiently well that the ultimate best-performing kernels are invariant across GA runs. This is the first justification to accept GA outputs as the “true” best-fitting kernels. Similarly, while the BIC metric is determined by the precise dataset, it ought to be sufficiently stable when the dataset undergoes a small modification, so that the top kernels can be interpreted as being the right ones for the population in question, and not just for the particular data subset picked.

The above robustness checks are summarized below and in Table 9. These confirm that the GA is stable both across its own runs (see the “Re-run” listings) and when subjected to slightly modifying the training dataset or “rolling” it in time. We observe that we recover essentially the same kernels (both in ${\mathcal{K}}_r$ and $\mathcal{K}_f$ ), and moreover the best kernels/hyperparameters are highly stable as we enlarge the dataset (see the “Robust” listings). This pertains both to the very top kernels, explicitly listed below, but also to the larger set of top-100 kernels, see the summary statistics in Table 9. Re-assuringly, the kernel lengthscales listed below change minimally when run on a larger dataset, confirming the stability of the MLE GP sub-routines. In particular, the same kernel is identified as the best one during the re-run in ${\mathcal{K}}_r$ , and it shows up yet again as the best for $\mathcal{D}_{\textrm{rob,1}}$ , with just slightly modified parameters:

\begin{align*}\textrm{original $\mathcal{D}$: } &0.4651 \cdot \textrm{M}52_a({37.7}) \cdot \textrm{M}52_y({52.2})\cdot \textrm{M}12_y({1821}) \cdot \textrm{M}12_c({7412}); \\\textrm{enlarged $\mathcal{D}_{\textrm{rob,1}}$: } &0.4646\cdot \textrm{M}52_a({37.7})\cdot \textrm{M}52_y({52.2})\cdot \textrm{M}12_y({1819})\cdot \textrm{M}12_c({7403}).\end{align*}

Four out of the five best kernels repeat when working with $\mathcal{D}_{\textrm{rob,1}}$ . This stability can be contrasted with Cairns et al. (Reference Cairns, Blake, Dowd, Coughlan, Epstein and Khalaf-Allah2011) who comment on sensitivity of SVD-based fitting to date range. The other alternatives continue to follow familiar substitution patterns. Of note, with the run over $\mathcal{D}_{\textrm{rob,1}}$ there is the appearance of $\textrm{Chy}_a$ , but no appearance of $\textrm{Min}_y, \textrm{M}32_y,$ or $\textrm{Meh}_c$ . As can be seen in Figure 2, $\textrm{Chy}_a$ is actually quite common.

Furthermore, all runs (original, re-run, enlarged dataset) always select $\textrm{M}52$ or Cauchy kernel for the Age effect, Matérn-1/2 (or sometimes Min) in Period, typically augmented with a smoother kernel like $\textrm{M}32_y,\textrm{Meh}_y,\textrm{Chy}_y,$ and $\textrm{M}12_c$ in Cohort. We furthermore record very similar frequency of different kernels among top-100 proposals, and similar BIC scores for the re-run.

Figure 5. Summary statistics of best kernels proposed by GA as a function of generation g.

As another validation of GA convergence, Figure 5 shows the evolution of fitness scores over generations. We display the BIC of the best kernel in generation g, as well as the second best (99% quantile across 200 kernels), 5th best (97.5%), 10th best (95%), and 20th best (90%) across the main run of the GA and a “re-run”. The experiments in each panel differ only through a different initial seed for the first generation. In both settings, only minimal performance increases (according to the minimum BIC) are found beyond generation 12 or so. Since in each new generation there is inherent randomness in newly proposed kernels, there is only distributional convergence of the BIC scores as new kernels are continuously tried out. This churn is indicated by the flat curves of the respective within-generation BIC quantiles. In sum, the GA converges to its “equilibrium” after about a dozen generations, validating our use of $G=20$ for analysis.

5.1.1. Robustness to prior mean

To examine the stability of top-ranking kernel compositions under different prior mean functions, we compare our primary choice of the linear prior mean $m_{lin}(\textbf{x}) = \beta_0 + \beta_{ag} x_{ag}$ to the following two alternatives:

1. (i) A constant mean function $m_c(\textbf{x}) = \beta_0$ ;

2. (ii) A calendar year-averaged mean $m_{ya}(\textbf{x}) = n_{\textrm{yr}}^{-1} \sum_{x_{\textrm{yr}} \in \mathcal{D}} y( x_{\textrm{ag}},x_{\textrm{yr}} )$ .

The constant mean $m_c$ is the simplest possible prior and drops the second coefficient $\beta_{ag}= 0$ , while still estimating $\beta_0$ alongside covariance hyperparameters during MLE. It is expected to yield a worse goodness-of-fit and also lead to longer lengthscales, as the GP must match the linear trend directly. The year-averaged mean $m_{ya}$ is an example of a pre-determined (rather than fitted) mean function that provides an age-specific data-driven trend; this would generally be a more accurate de-trending compared to the linear function in $m_{lin}({\cdot})$ . Year-averaging de-trending is used in Cairns et al. (Reference Cairns, Blake, Dowd, Coughlan, Epstein, Ong and Balevich2009) (comparing the Lee–Carter with and without cohort models), and also appears in the Renshaw–Haberman model (Renshaw and Haberman, Reference Renshaw and Haberman2006). Below our goal is not to compare performance across different mean functions but to assess robustness of covariance structure to choice of mean function.

Re-running the GA with these alternatives, we find the following top kernels:

\begin{align*}m_{ya}: \quad&0.1567\cdot \textrm{Chy}_a(45.4)\cdot (\textrm{RBF}_y(24.7) \textrm{M}12_y(579.9))\cdot \textrm{M}12_c(2562); \\m_c: \quad &2.1117\cdot \textrm{M}52_a(40.9)\cdot \textrm{M}12_y(7873)\cdot (\textrm{M}12_c(31494) \textrm{Meh}_c(0.473)).\end{align*}

In all, these results are remarkably similar to the primary run and demonstrate that the choice of the mean function mostly affects kernel lengthscales but not the selected kernel types. Throughout the three choices of $m({\cdot}),$ we consistently get a purely multiplicative kernel with four terms, including a smooth Age effect (captured either with Chy or $\textrm{M}52$ kernels that appear fully substitutable), a quasi-nonstationary Year effect with a $\textrm{M}12_y$ term, a similar cohort effect with a $\textrm{M}12_c$ term, and a fourth smooth term, mostly in Cohort using one of $\textrm{Meh},\textrm{Chy},\textrm{M}52$ kernels. We note that with a constant mean there is less of a mean-reversion effect, that is, the GP is more focused on extrapolating the Year pattern, see the very large lengthscales in Year and Cohort. Also, as expected, we find a reduced GP process variance $\sigma^2_{f,ya} \approx 0.15 \ll \sigma_{f,lin}^2 \approx 0.45$ for the year-averaged case (where the de-trended residuals modeled by the GP are smaller in magnitude) and a much larger $\sigma_{f,c}^2 \approx 2$ for the constant prior mean.

5.2. Analysis of model residuals

In Figure 6, the left panel displays the residuals that compare the realized log-mortality rates of JPN Females with the GP prediction from the best-performing kernel. The absence of any identifiable structure, especially along the SW-NE diagonals that correspond to Birth Cohorts, indicates a statistically sound fit, consistent with the expected uncorrelated and identically distributed residuals. Additionally, we observe distinct heteroskedasticity, where residuals for smaller Ages exhibit higher variance. This is due to the smaller number of deaths at those Ages, resulting in a more uncertain inferred mortality rate, despite the larger number of exposures. Generally, the observation variance is lowest around Age 80.

Figure 6. Left: residuals from the best kernel in $\mathcal{K}_f$ for JPN Females. Right: implied prior correlation $k(\textbf{x}_0, \textbf{x}^{\prime})$ of the best kernel as function of $\textbf{x}^{\prime}$ relative to the cell $\textbf{x}_0=(70,2010)$ shown as the red dot.

The right panel of Figure 6 shows the implied prior correlation relative to the cell (70, 2010). The strong diagonal shape indicates the importance of the cohort effect. Moreover, we observe that the correlation decays about the same in Period (vertical) as in Age (horizontal), with the inferred lengthscales imposing a dependence of about $\pm 12$ years in each direction.

5.3. Male versus female populations

We proceed to apply the GA to the JPN Male population. The rationale behind this comparison is the assumption of similar correlation structures between the genders, which enables us to both highlight the similarities and pinpoint the observed differences.

As expected, the JPN Male results (see Tables 11 and 12) strongly resemble those of JPN Females. Once again, we detect a strong indication of a single multiplicative term, characterized by APC structure with smooth Age and rough Period and Cohort effects. Just like for Females, a (rough) Cohort term is selected in all (100%) of the top-performing kernels for JPN Males. The best-fitting individual kernels, as shown in Table 12, are also similar, with M52 in Age, M52 or RBF in Year, and M12 in Cohort being identified as the optimal choices.

Table 11. Results from GA runs on JPN Male, US Male, and SWE Female. Throughout we search within the full set $\mathcal{K}_f$ . See Table 10 for the full definition of all the columns.

Table 12. Best-performing kernel in ${\mathcal{K}}_r$ and $\mathcal{K}_f$ for each of the four populations considered. $N_{pl}$ is the number of alternate kernels that have a BIC within 6.802 of the top kernel and hence are judged “plausible” based on the BF criterion.

Some differences, such as more Cohort-linked kernels for JPN Males versus Females, are also observed. The Cohort lengthscale is much larger for females (7600 vs. 2500), and so is the rough $\textrm{M}12_y$ lengthscale (1800 vs. 1100), while the Age lengthscales are almost identical. One interpretation is that there is more idiosyncratic noise in Male mortality, leading to faster correlation decay.

5.4. Analysis across countries

To offer a broader cross-section of the global mortality experience, we next also consider US males and Sweden females. In total, we thus analyze four datasets: US males, Japan females and males, and Sweden females. We note that Sweden is much smaller (10M population compared to 130M in Japan and 330M in USA) than the other two countries and therefore has much noisier data.

USA Males: The US male data lead to kernels of much higher length compared to all other datasets. The GA returns kernels with 5–6 base kernels and frequently includes two or even three additive terms. Moreover, the APC pattern is somewhat disrupted, possibly due to collinearity between the multiple Period and Cohort terms.

To demonstrate some of the observed characteristics, let us examine the top kernel in $\mathcal{K}_f$ , as presented in Table 12. This kernel comprises 11 terms, including 2 additive terms. However, we note that the second term has a significantly smaller coefficient, indicating that it serves as a “correction” term that is introduced to account for a less prominent and identifiable feature relative to the primary terms. Additionally, we observe that this kernel incorporates both rough Cohort term $\textrm{M}12_c$ and smoother ones, namely $\textrm{M}32_c$ and $\textrm{M}52_c$ . This points toward a multi-scale Cohort effect, where a few exceptional years (such as birth years during the Spanish Flu outbreak in 1918–1919) are combined with generational patterns (e.g., Baby Boomers vs. the Silent Generation). Unlike other datasets, the US Male data even include a $\textrm{RBF}_c$ term. Finally, the Age effect is described by the $\textrm{AR}2$ kernel, which is also commonly observed in the JPN Male population.

Moving down the list, there are also shorter kernels with a single component (no “+”), for example, $0.0129 \cdot \textrm{M}12_a(117.0) \cdot (\textrm{RBF}_y(18.3) \textrm{M}12_y(228.8)) \cdot (\textrm{M}32_c(27.1) \textrm{M}12_c(244.4))$ which is fourth-best, and the length-7 $0.0113 \cdot \textrm{M}12_a(88.7) \cdot (\textrm{M}32_y(20.3) \textrm{M}12_y(189.4))\cdot \textrm{M}12_c(183.3)$ which is seventh-best. In all, for US males we can find a plausible kernel of length 7, 9, 11, 13 when the best-performing one has length 11. This wide distribution of plausible kernel lengths (and a wide range of proposed kernel families) is illustrated in the right panel of Figure 3 and the middle column of Figure 7.

Figure 7. Frequency of appearance of different kernels from $\mathcal{K}_f$ in US, SWE, and JPN Male models.

Another sign that the US data have inherent complexity is the wide gap in BF of the best kernel in $\mathcal{K}_f$ compared to that in ${\mathcal{K}}_r$ , by far the biggest among all populations. Thus, restricting to ${\mathcal{K}}_r$ materially worsens the fit. In fact, we observe that all the top kernels in ${\mathcal{K}}_r$ are purely multiplicative (such as $\textrm{M}12_a \cdot \textrm{M}52_y \textrm{M}12_y \cdot \textrm{M}52_c \textrm{M}12_c$ ), which is unlikely to be the correct structure for these data and moreover hints at difficulty in capturing the correlation in each coordinate, leading to multiple Period and Cohort terms. Within $\mathcal{K}_f$ only 35 plausible alternatives are found, 3–5 times fewer than in other datasets.

SWE Female: The Swedish Females dataset turns out to have two distinguishing features. First, it yields the simplest and shortest kernels that directly match the APC structure of three multiplicative terms. The average kernel length reported in Table 11 is the smallest for SWE Females, and additive terms appear very rarely. When kernels with more than three terms are proposed, these are usually still all-multiplicative and add either a second Period term (e.g. $\textrm{Meh}_a \cdot \textrm{M}12_y \cdot \textrm{Meh}_y \cdot \textrm{Meh}_c$ ) or a second Age term, though both are smooth: $(\textrm{Chy}_a \textrm{RBF}_a) \cdot \textrm{M}12_y \cdot \textrm{Meh}_c$ ).

Second, and unlike all other populations above, the Cohort effect is ambiguous in Sweden. About 15% of the top performing kernels (and 30% in ${\mathcal{K}}_r$ ) have no Cohort terms at all, instead proposing two terms for the Period effect. Those that do include a Cohort effect, use either a $\textrm{RBF}_c$ or a $\textrm{Meh}_c$ term, indicating no short-term cohort features, but only generational ones. Nine of the top-10 kernels in $\mathcal{K}_f$ and only 7 out of 10 in ${\mathcal{K}}_r$ have Cohort terms. In contrast, in JPN and USA rough cohort terms are present in every single top-100 kernel. Once again, our results are consistent with the literature, see Murphy (Reference Murphy2010) who discusses the lack of clarity on cohort effect for the Swedish female population.

A third observation is that SWE Female shows a compression of BIC values, that is, a lot of different kernels are proposed with very similar BICs. The GA returns over two hundred kernels with a BF within a ratio of 30 to the top one. This indicates little evidence to distinguish many different choices from each other and could be driven by the lower complexity of the Swedish mortality data.

5.5. Discussion

Best kernel families: The barplots in Figure 7 suggest that there is no clear-cut covariance structure that fits mortality patterns. Consequently, the models often propose a combination (usually a product, sometimes a sum) of various kernels. Moreover, there is no one-size-fits-all solution as far as different populations are concerned. For instance, the Age effect is typically modeled via a $\textrm{M}12$ kernel for US Males, a $\textrm{M}32$ kernel for JPN Males, and a Chy (or RBF or $\textrm{M}52$ ) kernel for SWE Females. Additionally, Chy may be identified as a possible Period term for US and JPN Males but never for SWE Females. As such, it is recommended to select different kernels for different case studies. This represents one of the significant differences compared to the classical APC framework, where the SVD decomposition is invariant across datasets, and researchers must manually test numerous combinations, as illustrated in Cairns et al. (Reference Cairns, Blake, Dowd, Coughlan, Epstein and Khalaf-Allah2011).

Necessity of Cohort Effect: To assess the impact of including a Birth Cohort term, we re-run the GA while excluding all cohort-specific kernels. This is a straightforward adjustment to the implementation and can be used to test whether cohort effects could be adequately explained through a well-chosen Age–Period kernel combination.

We first evaluate our models by comparing the Bayesian information criterion (BIC) of the top kernel that excludes Cohort terms with that of the full $\mathcal{K}_f$ . Additionally, we examine the residuals heatmap to detect any discernible diagonal patterns. Our findings indicate that the Cohort effect is overwhelmingly needed for US Males, JPN Females, and JPN Males. The absence of a Cohort term results in a significant increase in BIC, with a difference of 235 for JPN Male, 198 for US Male, and 111 for JPN Female. To put things in perspective, a BIC difference is considered significant only when it surpasses 6.802. The results are confirmed by Figure 8, which illustrates pronounced diagonals in the bottom row (the no-cohort models), contradicting the assumption of independent residuals.

For Sweden, the difference is only 1.24, which corresponds to a BF of 0.2894, indicating no significance. Moreover, the corresponding no-cohort residuals still appear satisfactory, and the associated kernel $0.1125\cdot \textrm{Chy}_a(28.9) \cdot \textrm{Meh}_y(0.653)\cdot \textrm{M}12_y(1094)$ is ranked ninth-best in the original $\mathcal{K}_f$ . Once again, our results are consistent with the literature, see Murphy (Reference Murphy2010) who discusses the lack of clarity on cohort effect for Sweden.

Cairns et al. (Reference Cairns, Blake, Dowd, Coughlan, Epstein and Khalaf-Allah2011) suggested that cohort effects might be partially or completely explained by well-chosen age and period effects. We find a partial confirmation of this finding in that in many populations there are more than two kernels used to explain the Period and Cohort dependence, and there is a clear substitution between them. However, this is a second-order effect; the primary necessity of including Cohort is unambiguous except for Sweden. For Sweden females, the need for Birth Cohort dependence is quite weak.

Kernel Substitution: Through its mutation operations, the GA naturally highlights substitution effects among different kernel families. Substitution of one kernel with another is intrinsic to the evolution of the GA, and by ranking the kernels in terms of their BIC, we observe the presence of many compositions that achieve nearly same performance and differ just by one term. The above effect is especially noticeable in purely multiplicative kernels that are prevalent in all populations except the USA. In this case, we can frequently observe that one of the terms can be represented by two or three kernel families, with the rest of the terms staying fixed.

We observed that certain kernels are commonly used as substitutes for each other, such as the RBF and $\textrm{M}52$ kernels, as well as the $\textrm{M}12$ and Min kernels. Although the latter pair differs by stationarity, the sample paths generated by Min are visually indistinguishable from those generated by $\textrm{M}12$ when the lengthscale parameter $\ell$ is large (usually $\ell_{\textrm{Len}} > 1000$ ). These substitution patterns are in agreement with our synthetic results, as discussed in Section 4. As an example, when training on the Japan Female dataset and searching in ${\mathcal{K}}_r$ , the top kernel is of the form $\textrm{M}52_a \cdot (\textrm{RBF}_y \cdot \textrm{M}12_y) \cdot \textrm{M}12_c$ , while the second-best according to BIC is $\textrm{M}52_a \cdot (\textrm{M}52_y \cdot \textrm{M}12_y)\cdot \textrm{M}12_c$ . This preserves the same macro-structure while replacing one of the two Period terms, cf. Table 10. Additionally, the next two ranked kernels are very similar, but substitute $\textrm{M}12$ with Min: the third-best is $\textrm{M}52_a \cdot (\textrm{RBF}_y \cdot \textrm{Min}_y) \cdot \textrm{M}12_c$ , and the fourth-best is $\textrm{M}52_a \cdot (\textrm{M}52_y \cdot \textrm{Min}_y) \cdot \textrm{M}12_c$ . It is important to note that during substitution, the lengthscales (and sometimes process variances) change, as these have a different meaning for different kernel families. For example, the RBF lengthscales tend to be about 50% smaller than those for its substitute, $\textrm{M}52$ .

The Cauchy kernel and $\textrm{M}52$ are also interchangeable, although Chy and RBF are less so. This effect of multiple substitutes for smooth kernels is nicely illustrated for SWE Females, where the four top kernels fix the Period and Cohort effects according to $\textrm{M}12_y \cdot \textrm{Meh}_c$ and then propose any of $\textrm{Chy}_a, \textrm{RBF}_a, \textrm{Meh}_a, \textrm{M}52_a$ for the Age effect. Substitutions are more common within $\mathcal{K}_f$ , since the availability of more kernel families contributes to “collinearity” and hence more opportunities for substitution. At the same time, we occasionally observe the ability to find a more suitable kernel family in $\mathcal{K}_f$ . For example, often we observe both a rough and a smooth kernel in Year, indicating that neither M12 nor RBF fit well on their own; in that case $\textrm{M}32$ sometimes appears to be a better single substitute. Similarly, $\textrm{Meh}_c$ is the most common choice for SWE Females and is replaced with $\textrm{RBF}_c$ or $\textrm{RBF}_c \cdot \textrm{M}12_c$ in ${\mathcal{K}}_r$ . Consequently, proposed kernels from $\mathcal{K}_f$ tend to be a bit shorter on average than those from ${\mathcal{K}}_r$ .

A further substitution effect happens between Period and Cohort terms. In JPN Females and SWE Females, we tend to observe a total of three terms, and there is a substitution between using two Period and one Cohort or one Period and two Cohort terms. For instance, in JPN female among top-10 kernels we find both $\textrm{M}52_a \cdot \textrm{M}12_y \cdot (\textrm{Meh}_c \textrm{M}12_c)$ and $\textrm{M}52_a \cdot (\textrm{Meh}_y \textrm{M}12_y) \cdot \textrm{M}12_c$ .

Additivity and Nonstationarity: Returning to the topic of additive components, our results generally suggest that the additive structure is generally weak. Specifically, introducing an additional additive component often provides only a minor improvement in goodness of fit, which is offset by the complexity penalty in BIC, resulting in lower BICs for additive kernels. Hence, additive kernels tend to be rejected by the GA on the grounds of parsimony. As a result, most kernels with four terms (length 7) and the majority with five terms are purely multiplicative.

Summarizing nonstationarity is a challenging task due to various factors, so any table providing the percentage of nonstationary data should be viewed with reservation. A more comprehensive analysis of nonstationarity can be seen through the frequency diagrams shown in Figures 2 and 7. It is important to note that M12 lengthscales tend to be large in fitting. When M12 has a large lengthscale, the resulting processes are visually indistinguishable from those generated by the nonstationary Min kernel. Therefore, in our data, the presence of M12 indicates potential nonstationarity. From this perspective, our findings suggest that all populations exhibit a (potentially) nonstationary period effect. The other nonstationary kernels (Lin, Meh) are rare but do occur, mostly in the SWE female population.