Revisiting Whittaker-Henderson Smoothing

Origin

Whittaker-Henderson (WH) smoothing is a gradation method aimed at correcting the effect of sampling fluctuations on an observation vector. It is applied to evenly-spaced discrete observations. Initially proposed by Whittaker (1922) for constructing mortality tables and further developed by the works of Henderson (1924), it remains one of the most popular methods among actuaries for constructing experience tables in life insurance. Extending to two-dimensional tables, it can be used for studying various risks, including but not limited to: mortality, disability, long-term care, lapse, mortgage default, and unemployment.

The one-dimensional case

Let \(\mathbf{y}\) be a vector of observations and \(\mathbf{w}\) a vector of positive weights, both of size \(n\). The estimator associated with Whittaker-Henderson smoothing is given by:

\[ \hat{\mathbf{y}} = \underset{\boldsymbol{\theta}}{\text{argmin}}\{F(\mathbf{y},\mathbf{w},\boldsymbol{\theta}) + R_{\lambda,q}(\boldsymbol{\theta})\} \]

where:

• \(F(\mathbf{y},\mathbf{w},\boldsymbol{\theta}) = \underset{i = 1}{\overset{n}{\sum}} w_i(y_i - \theta_i)^2\) represents a fidelity criterion to the observations,

• \(R_{\lambda,q}(\boldsymbol{\theta}) = \lambda \underset{i = 1}{\overset{n - q}{\sum}} (\Delta^q\boldsymbol{\theta})_i^2\) represents a smoothness criterion.

In the latter expression, \(\Delta^q\) denotes the forward difference operator of order \(q\), such that for any \(i\in[1,n - q]\):

\[ (\Delta^q\boldsymbol{\theta})_i = \underset{k = 0}{\overset{q}{\sum}} \begin{pmatrix}q \\ k\end{pmatrix}(- 1)^{q - k} \theta_{i + k}. \]

Let us define \(W = \text{Diag}(\mathbf{w})\), the diagonal matrix of weights, and \(D_{n,q}\) as the order \(q\) difference matrix of dimensions \((n-q) \times n\), such that \((D_{n,q}\boldsymbol{\theta})_i = (\Delta^q\boldsymbol{\theta})_i\) for all \(i \in [1, n-q]\). The most commonly used difference matrices of order 1 and 2 have the following forms:

\[ D_{n,1} = \begin{bmatrix} -1 & 1 & 0 & \ldots & 0 \\ 0 & -1 & 1 & \ldots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \ldots & 0 & -1 & 1 \\ \end{bmatrix} \quad\text{and}\quad D_{n,2} = \begin{bmatrix} 1 & -2 & 1 & 0 & \ldots & 0 \\ 0 & 1 & -2 & 1 & \ldots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \ldots & 0 & 1 & -2 & 1 \\ \end{bmatrix}. \]

The fidelity and smoothness criteria can be rewritten in matrix form as:

\[ F(\mathbf{y},\mathbf{w},\boldsymbol{\theta}) = (\mathbf{y} - \boldsymbol{\theta})^TW(\mathbf{y} - \boldsymbol{\theta}) \quad \text{and} \quad R_{\lambda,q}(\boldsymbol{\theta}) = \lambda\boldsymbol{\theta}^TD_{n,q}^TD_{n,q}\boldsymbol{\theta}. \]

The associated estimator for smoothing becomes:

\[ \hat{\mathbf{y}} = \underset{\boldsymbol{\theta}}{\text{argmin}} \left\lbrace(\mathbf{y} - \boldsymbol{\theta})^TW(\mathbf{y} - \boldsymbol{\theta}) + \boldsymbol{\theta}^TP_{\lambda}\boldsymbol{\theta}\right\rbrace \]

where \(P_{\lambda} = \lambda D_{n,q}^TD_{n,q}\).

The two-dimensional case

In the bidimensional case, let us consider a matrix \(Y\) of observations and a matrix \(\Omega\) of non-negative weights, both of dimensions \(n_x \times n_z\). The estimator associated with the Whittaker-Henderson smoothing can be written as:

\[ \widehat{Y} = \underset{\Theta}{\text{argmin}}\{F(Y,\Omega, \Theta) + R_{\lambda,q}(\Theta)\} \]

where:

• \(F(Y,\Omega, \Theta) = \sum_{i = 1}^{n_x}\sum_{j = 1}^{n_z} \Omega_{i,j}(Y_{i,j} - \Theta_{i,j})^2\) represents a fidelity criterion to the observations,

• \(R_{\lambda,q}(\Theta) = \lambda_x \sum_{j = 1}^{n_z}\sum_{i = 1}^{n_x - q_x} (\Delta^{q_x}\Theta_{\bullet,j})_i^2 + \lambda_z \sum_{i = 1}^{n_x}\sum_{j = 1}^{n_z - q_z} (\Delta^{q_z}\Theta_{i,\bullet})_j^2\) is a smoothness criterion.

This latter criterion can be written as the sum of two one-dimensional regularization criteria, with orders \(q_x\) and \(q_z\), applied respectively to all rows and all columns of \(\Theta\). It is also possible to adopt matrix notations by defining \(\mathbf{y} = \textbf{vec}(Y)\), \(\mathbf{w} = \textbf{vec}(\Omega)\), and \(\boldsymbol{\theta} = \textbf{vec}(\Theta)\) as the vectors obtained by stacking the columns of the matrices \(Y\), \(\Omega\), and \(\Theta\), respectively. Additionally, let us denote \(W = \text{Diag}(\mathbf{w})\) and \(n = n_x \times n_z\). The fidelity and smoothness criteria can then be rewritten as linear combinations of the vectors \(\mathbf{y}\), \(\mathbf{w}\), and \(\boldsymbol{\theta}\):

\[ \begin{aligned} F(\mathbf{y},\mathbf{w}, \boldsymbol{\theta}) &= (\mathbf{y} - \boldsymbol{\theta})^TW(\mathbf{y} - \boldsymbol{\theta}) \\ R_{\lambda,q}(\boldsymbol{\theta}) &= \boldsymbol{\theta}^{T}(\lambda_x I_{n_z} \otimes D_{n_x,q_x}^{T}D_{n_x,q_x} + \lambda_z D_{n_z,q_z}^{T}D_{n_z,q_z} \otimes I_{n_x}) \boldsymbol{\theta}. \end{aligned} \]

and the estimator associated with the smoothing takes the same form as in the one-dimensional case:

\[ \hat{\mathbf{y}} = \underset{\boldsymbol{\theta}}{\text{argmin}} \left\lbrace(\mathbf{y} - \boldsymbol{\theta})^TW(\mathbf{y} - \boldsymbol{\theta}) + \boldsymbol{\theta}^TP_{\lambda}\boldsymbol{\theta}\right\rbrace \]

except in this case \(P_{\lambda} = \lambda_x I_{n_z} \otimes D_{n_x,q_x}^{T}D_{n_x,q_x} + \lambda_z D_{n_z,q_z}^{T}D_{n_z,q_z} \otimes I_{n_x}\).

Explicit solution

Whittaker-Henderson estimator has an explicit solution:

\[\hat{\mathbf{y}} = (W + P_{\lambda})^{-1}W\mathbf{y}.\]

Indeed, as a minimum, \(\hat{\mathbf{y}}\) satisfies:

\[0 = \left.\frac{\partial}{\partial \boldsymbol{\theta}}\right|_{\hat{\mathbf{y}}}\left\lbrace(\mathbf{y} - \boldsymbol{\theta})^{T}W(\mathbf{y} - \boldsymbol{\theta}) + \boldsymbol{\theta}^{T}P_{\lambda}\boldsymbol{\theta}\right\rbrace = - 2 X^TW(y - \hat{\mathbf{y}}) +2X^TP_{\lambda} \hat{\mathbf{y}}.\]

It follows that \((W + P_{\lambda})\hat{\boldsymbol{\theta}} = W\mathbf{y}\), and if \(W + P_{\lambda}\) is invertible, then \(\hat{\mathbf{y}}\) is indeed a solution for the original equation. When \(\lambda \neq 0\), \(W + P_{\lambda}\) is invertible as long as \(\mathbf{w}\) has \(q\) non-zero elements in the one-dimensional case, and \(\Omega\) has at least \(q_x \times q_z\) non-zero elements distributed over \(q_x\) different rows and \(q_z\) different columns in the two-dimensional case. These are sufficient conditions, which are always satisfied in practice and will not be demonstrated here.

How to obtain credibility intervals?

From the explicit solution of WH smoothing, \(\mathbb{E}(\hat{\mathbf{y}}) = (W + P_{\lambda})^{-1}W\mathbb{E}(\mathbf{y}) \ne \mathbb{E}(\mathbf{y})\) when \(\lambda \ne 0\). This implies that penalization introduces a smoothing bias, which prevents the construction of a confidence interval centered around \(\mathbb{E}(\mathbf{y})\). Therefore, in this section, we turn to a Bayesian approach where smoothing can be interpreted more naturally.

Let us suppose that \(\mathbf{y} | \boldsymbol{\theta}\sim \mathcal{N}(\boldsymbol{\theta}, W^{-})\) and \(\boldsymbol{\theta} \sim \mathcal{N}(0, P_{\lambda}^{-})\) where \(P_{\lambda}^{-}\) denotes the pseudo-inverse of the matrix \(P_{\lambda}\). The Bayes’ formula allows us to express the posterior likelihood \(f(\boldsymbol{\theta} | \mathbf{y})\) associated with these choices in the following form:

\[\begin{aligned} f(\boldsymbol{\theta} | \mathbf{y}) &= \frac{f(\mathbf{y} | \boldsymbol{\theta}) f(\boldsymbol{\theta})}{f(y)} \\ &\propto f(\mathbf{y} | \boldsymbol{\theta}) f(\boldsymbol{\theta}) \\ &\propto \exp\left(- \frac{1}{2}(\mathbf{y} - \boldsymbol{\theta})^{T}W(\mathbf{y} - \boldsymbol{\theta})\right)\exp\left(-\frac{1}{2}\boldsymbol{\theta}^TP_{\lambda} \boldsymbol{\theta}\right) \\ &\propto \exp\left(- \frac{1}{2}\left[(\mathbf{y} - \boldsymbol{\theta})^{T}W(\mathbf{y} - \boldsymbol{\theta}) + \boldsymbol{\theta}^TP_{\lambda}\boldsymbol{\theta}\right]\right). \end{aligned}\]

The mode of the posterior distribution, \(\hat{\boldsymbol{\theta}} = \underset{\boldsymbol{\theta}}{\text{argmax}} [f(\boldsymbol{\theta} | \mathbf{y})]\), also known as the maximum a posteriori (MAP) estimate, coincides with the explicit solution \(\hat{\mathbf{y}}\).

A second-order Taylor expansion of the log-posterior likelihood around \(\hat{\mathbf{y}} = \hat{\boldsymbol{\theta}}\) gives us:

\[\ln f(\boldsymbol{\theta} | \mathbf{y}) = \ln f(\hat{\boldsymbol{\theta}} | \mathbf{y}) + \left.\frac{\partial \ln f(\boldsymbol{\theta} | \mathbf{y})}{\partial \boldsymbol{\theta}}\right|_{\boldsymbol{\theta} = \hat{\boldsymbol{\theta}}}^{T}(\boldsymbol{\theta} - \hat{\boldsymbol{\theta}}) + \frac{1}{2}(\boldsymbol{\theta} - \hat{\boldsymbol{\theta}})^{T} \left.\frac{\partial^2 \ln f(\boldsymbol{\theta} | \mathbf{y})}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}^{T}}\right|_{\boldsymbol{\theta} = \hat{\boldsymbol{\theta}}}(\boldsymbol{\theta} - \hat{\boldsymbol{\theta}})\]

\[\text{where} \quad \left.\frac{\partial \ln f(\boldsymbol{\theta} | \mathbf{y})}{\partial \boldsymbol{\theta}}\right|_{\boldsymbol{\theta} = \hat{\boldsymbol{\theta}}} = 0 \quad \text{and} \quad \left.\frac{\partial^2 \ln f(\boldsymbol{\theta} | \mathbf{y})}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}^{T}}\right|_{\boldsymbol{\theta} = \hat{\boldsymbol{\theta}}} = - (W + P_{\lambda}).\]

Note that this last derivative no longer depends on \(\boldsymbol{\theta}\), and higher-order derivatives beyond 2 are all zero. The Taylor expansion allows for an exact computation of \(\ln f(\boldsymbol{\theta} | \mathbf{y})\). By substituting the result back into the Taylor expansion, we obtain:

\[\begin{aligned} f(\boldsymbol{\theta} | \mathbf{y}) &\propto \exp\left[\ln f(\hat{\boldsymbol{\theta}} | \mathbf{y}) - \frac{1}{2} (\boldsymbol{\theta} - \hat{\boldsymbol{\theta}})^{T}(W + P_{\lambda})(\boldsymbol{\theta} - \hat{\boldsymbol{\theta}})\right] \\ &\propto \exp\left[- \frac{1}{2} (\boldsymbol{\theta} - \hat{\boldsymbol{\theta}})^{T}(W + P_{\lambda})(\boldsymbol{\theta} - \hat{\boldsymbol{\theta}})\right] \end{aligned}\]

which can immediately be recognized as the density of the \(\mathcal{N}(\hat{\boldsymbol{\theta}},(W + P_{\lambda})^{- 1})\) distribution.

The assumption \(\boldsymbol{\theta} \sim \mathcal{N}(0, P_{\lambda}^{-})\) corresponds to a simple Bayesian formalization of the smoothness criterion. It reflects a (improper) prior belief of the modeler about the underlying distribution of the observation vector \(\mathbf{y}\).

The use of Whittaker-Henderson smoothing in the Bayesian framework and the construction of credible intervals are conditioned on the validity of the assumption \(\mathbf{y} | \boldsymbol{\theta}\sim \mathcal{N}(\boldsymbol{\theta}, W^{-})\). The components of the observation vector should be independent and have known variances. The weight vector \(\mathbf{w}\) used should contain the inverses of these variances. Under these assumptions, credible intervals at \(100(1 - \alpha)\)% for smoothing can be constructed and take the form:

\[\mathbb{E}(\mathbf{y}) | \mathbf{y} \in \left[\hat{\mathbf{y}} \pm \Phi\left(1 -\frac{\alpha}{2}\right)\sqrt{\textbf{diag}\left\lbrace(W + P_{\lambda})^{-1}\right\rbrace}\right]\]

with probability \(1 -\frac{\alpha}{2}\) where \(\hat{\mathbf{y}} = (W + P_{\lambda})^{- 1}W\mathbf{y}\) and \(\Phi\) denotes the cumulative distribution function of the standard normal distribution. According to Marra and Wood (2012), these credible intervals provide satisfactory coverage of the corresponding frequentist confidence intervals and can therefore be used as practical substitutes.

Duration models framework: one-dimensional case

Consider the observation of \(m\) individuals in a longitudinal study subject to left truncation and right censoring phenomena. Suppose we want to estimate a single distribution that depends on only one continuous explanatory variable, denoted by \(x\). For illustration purposes, we can consider it as a mortality distribution with the explanatory variable of interest \(x\) representing age. Such a distribution is fully characterized by one of the following quantities:

• the cumulative distribution function \(F(x)\) or its complement, the survival function \(S(x) = 1 - F(x)\),
• the associated probability density function \(f(x) = - \frac{\text{d}}{\text{d}x}S(x)\),
• the instantaneous hazard function \(\mu(x) = - \frac{\text{d}}{\text{d}x}\ln S(x)\).

Suppose that the considered distribution depends on a vector of parameters \(\boldsymbol{\theta}\) that we want to estimate using maximum likelihood. The likelihood associated with the observation of the individuals can be written as follows:

\[ \mathcal{L}(\boldsymbol{\theta}) = \underset{i = 1}{\overset{m}{\prod}} \left[\frac{f(x_i + t_i,\boldsymbol{\theta})}{S(x_i,\boldsymbol{\theta})}\right]^{\delta_i}\left[\frac{S(x_i + t_i,\boldsymbol{\theta})}{S(x_i,\boldsymbol{\theta})}\right]^{1 - \delta_i} \]

Where \(x_i\) represents the age at the start of observation, \(t_i\) represents the observation duration, i.e., the time elapsed between the start and end dates of observation, and \(\delta_i\) is the indicator of event observation, which takes the value 1 if the event of interest is observed and 0 if the observation is censored. We will not go into the details of calculating these three quantities, which should take into account the individual-specific information such as the subscription date, redemption date if applicable, as well as the global characteristics of the product, such as the presence of a deductible or medical selection phenomenon, and the choice of a restricted observation period due to data quality issues or delays in the reporting of observations. These factors typically lead to a narrower observation period than the actual presence period of individuals in the portfolio.

The various quantities introduced above are related by the following relationships:

\[ S(x) = \exp\left(\underset{u = 0}{\overset{x}{\int}}\mu(u)\text{d}u\right)\quad \text{and} \quad f(x) = \mu(x)S(x). \]

The maximization of the likelihood in that equation is equivalent to the maximization of the associated log-likelihood, which can be rewritten using only the instantaneous hazard function (also known as the mortality rate in the case of the death risk):

\[ \ell(\boldsymbol{\theta}) = \underset{i = 1}{\overset{m}{\sum}} \left[\delta_i \ln\mu(x_i + t_i,\boldsymbol{\theta}) - \underset{u = 0}{\overset{t_i}{\int}}\mu(x_i + u,\boldsymbol{\theta})\text{d}u\right] \]

To discretize the problem, let us assume that the mortality rate is piecewise constant over one-year intervals between two integer ages. Formally, we have \(\mu(x + \epsilon) = \mu(x)\) for all \(x \in \mathbb{N}\) and \(\epsilon \in [0,1[\). Furthermore, if \(\mathbf{1}\) denotes the indicator function, then for any \(0 \leq a < x_{\max}\), we have \(\sum_{x = x_{\min}}^{x_{\max}} \mathbf{1}(x \leq a < x + 1) = 1\), where \(x_{\min} = \min(\mathbf{x})\) and \(x_{\max} = \max(\mathbf{x})\). The log-likelihood can be rewritten as:

\[ \begin{aligned} \ell(\boldsymbol{\theta}) = \underset{i = 1}{\overset{m}{\sum}} &\left[\underset{x = x_{\min}}{\overset{x_{\max}}{\sum}} \delta_i\mathbf{1}(x \le x_i + t_i < x + 1)\ln\mu(x_i + t_i,\boldsymbol{\theta})\right. \\ &- \left.\underset{u = 0}{\overset{t_i}{\int}}\underset{x = x_{\min}}{\overset{x_{\max}}{\sum}} \mathbf{1}(x \le x_i + u < x + 1)\mu(x_i + u,\boldsymbol{\theta})\text{d}u\right]. \end{aligned} \]

The assumption of piecewise constant mortality rate implies that:

\[ \begin{aligned} \mathbf{1}(x \le x_i + t_i < x + 1) \ln\mu(x_i + t_i,\boldsymbol{\theta}) &= \mathbf{1}(x \le x_i + t_i < x + 1) \ln\mu(x,\boldsymbol{\theta})\quad\text{and}\\ \mathbf{1}(x \le x_i + u < x + 1)\mu(x_i + u,\boldsymbol{\theta}) &= \mathbf{1}(x \le x_i + u < x + 1) \ln\mu(x,\boldsymbol{\theta}). \end{aligned} \]

It is then possible to interchange the two summations to obtain the following expressions:

\[ \begin{aligned} \ell(\boldsymbol{\theta}) &= \underset{x = x_{\min}}{\overset{x_{\max}}{\sum}} \left[\ln\mu(x,\boldsymbol{\theta}) d(x) - \mu(x,\boldsymbol{\theta}) e_c(x)\right] \quad \text{where} \\ d(x) & = \underset{i = 1}{\overset{m}{\sum}} \delta_i \mathbf{1}(x \le x_i + t_i < x + 1) \quad \text{and} \\ e_c(x) & = \underset{i = 1}{\overset{m}{\sum}}\underset{u = 0}{\overset{t_i}{\int}}\mathbf{1}(x \le x_i + u < x + 1)\text{d}u = \underset{i = 1}{\overset{m}{\sum}} \left[\min(t_i, x - x_i + 1) - \max(0, x - x_i)\right]^+ \end{aligned} \]

by denoting \(a^+ = \max(a, 0)\), where \(d(x)\) and \(e_c(x)\) correspond to the number of observed deaths between ages \(x\) and \(x + 1\) and the sum of observation durations of individuals between these ages, respectively (the latter quantity is also known as central exposure to risk).