* Section 4 introduced generalized WH smoothing, in which the penalization is applied directly to the likelihood rather than a normal approximation. This refined method yields more accurate results, especially in situations where the available data volume is limited but the number of combinations is high, such as in the two-dimensional case.

* Section 5 advocated for smoothing parameter selection via marginal likelihood (or its Laplace approximation, LAML), offering a principled and robust alternative to heuristic criteria like AIC or GCV.

* Section 6 presented two computational improvements: one exploits the banded structure of WH matrices to reduce runtime by up to a factor of 25; the other relies on reduced-rank smoothing basis, leading to even faster estimation (up to 250 × speed-up) with limited loss in accuracy, slightly outperforming P-splines.

* Section 7 addressed extrapolation: while WH smoothing naturally extends beyond the data range, constraints are needed in two dimensions to preserve the original fit. We proposed a method to extrapolate while accounting for both structural uncertainty and innovation error, and provided credible intervals accordingly.

All these techniques are available in the WH package for the statistical software R (R Core Team 2025), including automated smoothing parameter selection and constrained extrapolation with uncertainty quantification.

## Limitations and outlook

Despite its strong practical appeal, WH smoothing has limitations that suggest several avenues for future work:

* Regular spacing requirement: WH smoothing assumes evenly spaced observations, which aligns well with standard life insurance grids (age and/or duration). However, this is less suitable when events are concentrated in a short period, such as in disability or long-term care claims. One solution is to combine finer discretization in early durations with methods like P-splines that accommodate irregular grids. Alternatively, and adaptive WH smoothing procedure (based on the ideas in Ruppert and Carroll 2000; Krivobokova, Crainiceanu, and Kauermann 2008) could offer a way to retain regular spacing while varying the smoothness locally.

* Limited covariate handling: The basic WH framework does not accommodate additional explanatory variables (e.g., gender or policy features). However, WH smoothing can be extended using ideas from smoothing spline ANOVA and hierarchical models (Lee and Durban 2011; Gu 2013), allowing for structured random effects and flexible interactions. This opens the door to richer, more personalized experience modelling while preserving interpretability.
In sum, revisiting WH smoothing through a modern lens reinforces its theoretical foundations and offers practitioners fast, transparent, and adaptable tools for experience modelling. It remains a compelling alternative to more recent—yet often more opaque—techniques when working with evenly spaced discrete data.

## References

Akaike, Hirotsugu. 1973. "Information Theory and an Extension of the Maximum Likelihood Principle." In *2nd International Symposium on Information Theory, 1973*.

Anderssen, RS, and Peter Bloomfield. 1974. "A Time Series Approach to Numerical Differentiation." *Technometrics* 16 (1): 69-75.

Biessy, Guillaume. 2015. "Long-Term Care Insurance: A Multi-State Semi-Markov Model to Describe the Dependency Process in Elderly People." *Bulletin Français d'Actuariat* 15(29): 41-73.

———. 2016. "Semi-Markov Modeling of the Loss of Autonomy Among Elderly People: Application to Long-Term Care Insurance." PhD thesis, Paris Saclay.

Bohlmann, Georg. 1899. "Ein Ausgleichungsproblem." *Nachrichten von Der Gesellschaft Der Wissenschaften Zu Göttingen, Mathematisch-Physikalische Klasse* 1899: 260-71.

Brent, Richard P. 1973. "Algorithms for Minimization Without Derivatives, Chap. 4." Prentice-Hall, Englewood Cliffs, NJ.

Brooks, RJ, M Stone, FY Chan, and LK Chan. 1988. "Cross-Validatory Graduation." *Insurance: Mathematics and Economics* 7 (1): 59-66.

Canadian Institute of Actuaries. 2017. "Final Report: Canadian Mortality Improvement Model - MI-2017." Canadian Institute of Actuaries. https://www.cia-ica.ca/app/themes/wicket/custom/dl_file.php?fid=14926&p=35281.

Carballo, Alba, Maria Durban, Göran Kauermann, and Dae-Jin Lee. 2021. "A General Framework for Prediction in Penalized Regression." *Statistical Modelling* 21 (4): 293-312.

Carballo, Alba, Maria Durban, and Dae-Jin Lee. 2021. "Out-of-Sample Prediction in Multidimensional p-Spline Models." *Mathematics* 9 (15): 1761.

Clifford, Peter. 1977. "Nonidentifiability in Stochastic Models of Illness and Death." *Proceedings of the National Academy of Sciences* 74 (4): 1338-40.

Cornea-Madeira, Adriana. 2017. "The Explicit Formula for the Hodrick-Prescott Filter in a Finite Sample." *Review of Economics and Statistics* 99 (2): 314-18.

Currie, Iain D. 2013. "Smoothing Constrained Generalized Linear Models with an Application to the Lee-Carter Model." *Statistical Modelling* 13 (1): 69-93.

Currie, Iain D, Maria Durban, and Paul HC Eilers. 2004. "Smoothing and Forecasting Mortality Rates." *Statistical Modelling* 4 (4): 279-98.

———. 2006. "Generalized Linear Array Models with Applications to Multidimensional Smoothing." *Journal of the Royal Statistical Society: Series B (Statistical Methodology)* 68 (2): 259-80.
Delwarde, Antoine, Michel Denuit, and Paul Eilers. 2007. "Smoothing the Lee-Carter and Poisson Log-Bilinear Models for Mortality Forecasting: A Penalized Log-Likelihood Approach." *Statistical Modelling* 7 (1): 29-48.

Demmler, A, and Christian Reinsch. 1975. "Oscillation Matrices with Spline Smoothing." *Numerische Mathematik* 24 (5): 375-82.

Dempster, Arthur P, Nan M Laird, and Donald B Rubin. 1977. "Maximum Likelihood from Incomplete Data via the EM Algorithm." *Journal of the Royal Statistical Society: Series B (Methodological)* 39 (1): 1-22.

Eilers, Paul H. C., and Brian D. Marx. 1996. "Flexible Smoothing with B-Splines and Penalties." *Statistical Science* 11 (2): 89-102.

Fellner, William H. 1986. "Robust Estimation of Variance Components." *Technometrics* 28 (1): 51-60.

Fix, Evelyn, and Jerzy Neyman. 1951. "A Simple Stochastic Model of Recovery, Relapse, Death and Loss of Patients." *Human Biology* 23 (3): 205-41.

Giesecke, Lee, and Defense Manpower Data Center. 1981. "Use of the Chi-Square Statistic to Set Whittaker-Henderson Smoothing Coefficients." Smoothing.

Golub, Gene H, and Charles F Van Loan. 2013. *Matrix Computations*. JHU press.

Gschlössl, Susanne, Pascal Schoenmaekers, and Michel Denuit. 2011. "Risk Classification in Life Insurance: Methodology and Case Study." *European Actuarial Journal* 1: 23-41.

Gu, Chong. 1992. "Cross-Validating Non-Gaussian Data." *Journal of Computational and Graphical Statistics* 1 (2): 169-79.

———. 2013. *Smoothing Spline ANOVA Models*. Vol. 297. Springer.

Hastie, Trevor J, and Robert J Tibshirani. 1990. *Generalized Additive Models*. Vol. 43. CRC press.

Henderson, Robert. 1924. "A New Method of Graduation." *Transactions of the Actuarial Society of America* 25: 29-40.

Hodrick, Robert J, and Edward C Prescott. 1997. "Postwar US Business Cycles: An Empirical Investigation." *Journal of Money, Credit, and Banking*, 1-16.

Hoem, Jan M. 1971. "Point Estimation of Forces of Transition in Demographic Models." *Journal of the Royal Statistical Society: Series B (Methodological)* 33 (2): 275-89.

Kauermann, Göran. 2005. "A Note on Smoothing Parameter Selection for Penalized Spline Smoothing." *Journal of Statistical Planning and Inference* 127 (1-2): 53-69.

Knorr, Frank E. 1984. "Multidimensional Whittaker-Henderson Graduation." *Transactions of Society of Actuaries* 36: 213-55.

Krivobokova, Tatyana, Ciprian M Crainiceanu, and Göran Kauermann. 2008. "Fast Adaptive Penalized Splines." *Journal of Computational and Graphical Statistics* 17 (1): 1-20.

Lee, Dae-Jin, and Maria Durban. 2011. "P-Spline ANOVA-Type Interaction Models for Spatio-Temporal Smoothing." *Statistical Modelling* 11 (1): 49-69.

Marra, Giampiero, and Simon N Wood. 2012. "Coverage Properties of Confidence Intervals for Generalized Additive Model Components." *Scandinavian Journal of Statistics* 39 (1): 53-74.

Nelder, John Ashworth, and Roger Mead. 1965. "A Simplex Method for Function Minimization." *The Computer Journal* 7 (4): 308-13.
Nelder, John Ashworth, and Robert WM Wedderburn. 1972. "Generalized Linear Models." *Journal of the Royal Statistical Society: Series A (General)* 135 (3): 370-84.

Patterson, H. D., and R. Thompson. 1971. "Recovery of Inter-Block Information When Block Sizes Are Unequal." *Biometrika* 58: 545-54.

R Core Team. 2025. *R: A Language and Environment for Statistical Computing*. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Reinsch, Christian H. 1967. "Smoothing by Spline Functions." *Numerische Mathematik* 10 (3): 177-83.

Reiss, Philip T, and R Todd Ogden. 2009. "Smoothing Parameter Selection for a Class of Semiparametric Linear Models." *Journal of the Royal Statistical Society: Series B (Statistical Methodology)* 71 (2): 505-23.

Rodriguez-Alvarez, Maria Xosé, Dae-Jin Lee, Thomas Kneib, Maria Durban, and Paul Eilers. 2015. "Fast Smoothing Parameter Separation in Multidimensional Generalized p-Splines: The SAP Algorithm." *Statistics and Computing* 25 (5): 941-57.

Rodríguez-Álvarez, María Xosé, Maria Durban, Dae-Jin Lee, and Paul HC Eilers. 2019. "On the Estimation of Variance Parameters in Non-Standard Generalised Linear Mixed Models: Application to Penalised Smoothing." *Statistics and Computing* 29: 483-500.

Ruppert, David, and Raymond J Carroll. 2000. "Theory & Methods: Spatially-Adaptive Penalties for Spline Fitting." *Australian & New Zealand Journal of Statistics* 42 (2): 205-23.

Schall, Robert. 1991. "Estimation in Generalized Linear Models with Random Effects." *Biometrika* 78 (4): 719-27.

Society of Actuaries. 2018. "A Practitioner's Guide to Statistical Mortality Graduation." https://www.soa.org/globalassets/assets/files/resources/tables-calcs-tools/2018-stat-mort-graduation.pdf.

Taylor, Greg. 1992. "A Bayesian Interpretation of Whittaker-Henderson Graduation." *Insurance: Mathematics and Economics* 11 (1): 7-16.

Verrall, R.J. 1993. "A State Space Formulation of Whittaker Graduation, with Extensions." *Insurance: Mathematics and Economics* 13 (1): 7-14.

Wahba, Grace. 1980. *Spline Bases, Regularization, and Generalized Cross Validation for Solving Approximation Problems with Large Quantities of Noisy Data*. University of Wisconsin.

———. 1985. "A Comparison of GCV and GML for Choosing the Smoothing Parameter in the Generalized Spline Smoothing Problem." *The Annals of Statistics*, 1378-1402.

Weinert, Howard L. 2007. "Efficient Computation for Whittaker-Henderson Smoothing." *Computational Statistics & Data Analysis* 52 (2): 959-74.

Whittaker, Edmund Taylor. 1923. "On a New Method of Graduation." *Proceedings of the Edinburgh Mathematical Society* 41: 63-75.

Wood, Simon N. 2011. "Fast Stable Restricted Maximum Likelihood and Marginal Likelihood Estimation of Semiparametric Generalized Linear Models." *Journal of the Royal Statistical Society: Series B (Statistical Methodology)* 73 (1): 3-36.

———. 2017. *Generalized Additive Models: An Introduction with r*. 2nd ed. chapman; hall/CRC.

Wood, Simon N, and Matteo Fasiolo. 2017. "A Generalized Fellner-Schall Method for Smoothing Parameter Optimization with Application to Tweedie Location, Scale and Shape Models."
Biometrics 73 (4): 1071–81.

Wood, Simon N, Zheyuan Li, Gavin Shaddick, and Nicole H Augustin. 2017. “Generalized Additive Models for Gigadata: Modeling the UK Black Smoke Network Daily Data.” Journal of the American Statistical Association 112 (519): 1199–1210.

# Appendices

## A Exposure computation in the survival analysis framework

This appendix outlines the methodology used to compute central exposure to risk in survival analysis, for both univariate and bivariate settings.

The derivation relies on standard assumptions of left truncation and non-informative right censoring, and models the hazard rate as piecewise constant over one-year intervals. This discretization allows reformulating the continuous-time log-likelihood in terms of aggregated death counts and exposure durations.

In the one-dimensional case, exposure corresponds to the total time under observation within each integer age interval. The same principle extends naturally to the two-dimensional case, such as mortality modelling by age and duration in LTC, where exposure is computed over the age-duration grid using the same piecewise constant assumption.

These quantities—death counts and central exposures—form the inputs to the generalized Whittaker-Henderson smoothing approach used throughout the paper.

### One-dimensional case

Consider the observation of *m* individuals in a longitudinal study subject to left truncation and non-informative right censoring. Suppose we aim to estimate a distribution that depends on only one continuous explanatory variable, denoted by *x*. One may for example think of a mortality distribution with the explanatory variable of interest *x* representing age. Such a distribution is fully characterized by either 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{d}{dx}S(x)$,

* the instantaneous hazard function $\mu(x) = -\frac{d}{dx}\ln S(x)$.

Those 3 quantities are related by the following relationships:

$$S(x) = \exp \left( \int_{u=0}^{x} \mu(u) du \right) \quad \text{and} \quad f(x) = \mu(x)S(x). \qquad (15)$$
Suppose that the considered distribution depends on a vector of parameters $\theta$ estimated using maximum likelihood. The likelihood associated with the observation of the individuals takes the form:

$$ \mathcal{L}(\theta) = \prod_{i=1}^{m} \left[ \frac{f(x_i + t_i, \theta)}{S(x_i, \theta)} \right]^{\delta_i} \left[ \frac{S(x_i + t_i, \theta)}{S(x_i, \theta)} \right]^{1-\delta_i} \quad (16) $$

where $x_i$ represents the age at the start of observation, $t_i$ represents the duration of observation for individual $i$ 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 instead censored. We will not go into the details of how these three quantities are derived, however they should take into account individual-specific information such as the subscription date, lapse date if applicable, as well as the global characteristics of the product such as the presence of a waiting period or medical selection phenomenon, and the choice of a restricted observation period due to delays in the reporting of event of interests. These factors typically result in a narrower effective observation window than the actual time individuals spend in the portfolio.

Using Equation 15, the log-likelihood associated with Equation 16 can be rewritten using only the instantaneous hazard function (also known as force of mortality in the case of the death risk):

$$ \ell(\theta) = \sum_{i=1}^{m} \left[ \delta_i \ln \mu(x_i + t_i, \theta) - \int_{u=0}^{t_i} \mu(x_i + u, \theta) du \right] \quad (17) $$

We discretize the problem by assuming that the mortality rate is piecewise constant over one-year intervals between two integer ages or more formally $\mu(x + \epsilon) = \mu(x)$ for all $x \in \mathbb{N}$ and $\epsilon \in [0, 1[$. Further note that, if $\mathbf{1}$ denotes the indicator function, then for any $x_{\min} \le a < x_{\max}$, we have $\sum_{x=x_{\min}}^{x_{\max}} \mathbf{1}(x \le a < x+1) = 1$, where $x_{\min} = \min(\mathbf{x})$ and $x_{\max} = \max(\mathbf{x})$. Equation 17 may therefore be rewritten as:

$$ \ell(\theta) = \sum_{i=1}^{m} \left[ \sum_{x=x_{\min}}^{x_{\max}} \delta_i \mathbf{1}(x \le x_i + t_i < x+1) \ln \mu(x_i + t_i, \theta) - \int_{u=0}^{t_i} \sum_{x=x_{\min}}^{x_{\max}} \mathbf{1}(x \le x_i + u < x+1) \mu(x_i + u, \theta) du \right]. $$

The assumption of piecewise constant mortality rates implies that:

$$ \mathbf{1}(x \le x_i + t_i < x+1) \ln \mu(x_i + t_i, \theta) = \mathbf{1}(x \le x_i + t_i < x+1) \ln \mu(x, \theta) \quad \text{and} \\ \mathbf{1}(x \le x_i + u < x+1) \mu(x_i + u, \theta) = \mathbf{1}(x \le x_i + u < x+1) \mu(x, \theta). $$

It is then possible to interchange the two summations to obtain the following expressions:
$$ \ell(\boldsymbol{\theta}) = \sum_{x=x_{\min}}^{x_{\max}} [\ln \mu(x, \boldsymbol{\theta})d(x) - \mu(x, \boldsymbol{\theta})e_c(x)] \quad \text{where} $$

$$ d(x) = \sum_{i=1}^{m} \delta_i \mathbf{1}(x \le x_i + t_i < x + 1) \quad \text{and} $$

$$ e_c(x) = \sum_{i=1}^{m} \int_{u=0}^{t_i} \mathbf{1}(x \le x_i + u < x + 1) du = \sum_{i=1}^{m} [\min(t_i, x - x_i + 1) - \max(0, x - x_i)]^+ $$

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.

## Two-dimensional case

The extension of the proposed approach to the two-dimensional framework requires only minor adjustments to the previous reasoning. Let $z_{\min} = \min(\mathbf{z})$ and $z_{\max} = \max(\mathbf{z})$. The piecewise constant assumption for the mortality rate needs to be extended to the second dimension. Formally, assume that $\mu(x + \epsilon, z + \xi) = \mu(x, z)$ for all pairs $x, z \in \mathbb{N}$ and $\epsilon, \xi \in [0, 1[$. The sums involving the variable $x$ are then replaced by double sums considering all combinations of $x$ and $z$. The log-likelihood becomes:

$$ \ell(\boldsymbol{\theta}) = \sum_{x=x_{\min}}^{x_{\max}} \sum_{z=z_{\min}}^{z_{\max}} [\ln \mu(x, z, \boldsymbol{\theta})d(x, z) - \mu(x, z, \boldsymbol{\theta})e_c(x, z)] \quad \text{where} $$

$$ d(x, z) = \sum_{i=1}^{m} \delta_i \mathbf{1}(x \le x_i + t_i < x + 1) \mathbf{1}(z \le z_i + t_i < z + 1) \quad \text{and} $$

$$ \mathbf{e}_c(x, z) = \sum_{i=1}^{m} \int_{u=0}^{t_i} \mathbf{1}(x \le x_i + u < x + 1) \mathbf{1}(z \le z_i + u < z + 1) du \\ = \sum_{i=1}^{m} [\min(t_i, x + 1 - x_i, z + 1 - z_i) - \max(0, x - x_i, z - z_i)]^+ . $$

## B Simulated datasets

This appendix details the simulation process used to generate the datasets on which the comparative analysis of accuracy and computational efficiency is based.
**General approach**

The simulated datasets were constructed following a five-step methodology:

1. Define hypothetical underlying laws to serve as the ground truth.

2. Generate an initial population of insured individuals.

3. Simulate life outcomes for each individual.

4. Extract samples of a predetermined size from the simulated population.

5. Compute aggregated event counts and central exposures for each sample.

**Defining the underlying laws**

The synthetic laws used in the simulations are not meant to be accurate representations of real-world phenomena. However, they incorporate features commonly observed in mortality and long-term care (LTC) experience studies to ensure plausible dynamics for testing purposes.

**General population mortality**

Mortality rates are derived from the Human Mortality Database for France in 2019. The dataset includes death counts and central exposure to risk by age (0 to 109) and gender. A smoothing algorithm is applied to reduce sampling noise, particularly at younger ages. The resulting rates define the general population mortality.

**Insured population mortality**

To reflect the well-documented observation that insured individuals generally exhibit lower mortality than the general population—especially at younger ages—we apply a smooth logistic adjustment factor to the general population mortality. This factor transitions from 40% at age 30 to 100% at age 90, with a midpoint of 70% at age 60. The adjusted rates define the insured population mortality, which is used for all annuity simulations.
## Long-Term Care (LTC) transition and mortality laws

For LTC simulations, assumptions are needed for:

* Autonomous mortality (i.e., mortality of individuals not in LTC),

* Incidence of entry into LTC,

* Mortality within LTC.

We draw upon assumptions from a technical note published by the reinsurer SCOR in 1995, adapted to use the previously defined insured population mortality $q_{\text{ref}}(x, g)$ for age $x$ and gender $g$:

* $q_a(x, g) = 0.8 \times q_{\text{ref}}(x, g)$ (autonomous mortality),

* $i(x) = 5.535 \times 10^{-3} \exp[(x - 52)/8]$ (LTC incidence),

* $q_i(x, g) = 2 \times q_{\text{ref}}(x, g) + 0.035$ (initial LTC mortality).

We refine the disabled mortality law to include a shock at LTC onset, reflecting heightened mortality during the first years in LTC, caused by cancer-related admissions (see for example Biessy 2015 for details about this phenomenon and impacts on curves for mortality in LTC). The refined mortality is:

$$q_i(x, t, g) = 2 \times q_{\text{ref}}(x, g) + 0.035 + K(g)f(x - t)h(t).$$

where:

* $t$ is the time since onset of LTC,

* $K(g)$ encodes a gender-specific intensity (0.5 for females, 0.75 for males),

* $f$ and $h$ are logistic-shaped modifiers to reflect attenuation with age and time since entry, respectively.

This formulation captures the elevated initial mortality due to severe conditions like terminal cancers, which tapers off within two years or by age 90. For further discussion, see Biessy (2016).
## Simulating a population of insured lives

We simulate 1,000,000 insured individuals over a subscription period from 1990 to 2009. Eligibility spans:

* ages 20 to 65 for annuity policies;

* ages 50 to 75 for LTC policies.

Subscription age, year, and gender are drawn with replacement from weighted distributions based on French population demographics. Birth dates and subscription dates are assigned uniformly at random within valid ranges.

## Simulating life trajectories

For each individual, we simulate life outcomes from their subscription date up to December 31, 2024. The process involves:

1. Dividing the time axis into intervals delimited by birthdays.

2. For each interval, computing event probabilities based on applicable mortality or incidence rates.

3. Drawing events using uniform random variables:

* In the LTC case, distinguishing between autonomous death and LTC incidence.

* Recording the event date accordingly.

4. If no terminating event occurs, proceeding to the next interval.

For individuals who enter LTC, a second simulation phase begins, spanning from LTC entry to the end of observation. This time, the timeline is segmented by both age and duration-in-LTC anniversaries. Within each subperiod, we compute and simulate death-in-LTC events.

## Sampling subsets from the simulated population

### Sample from the simulated population to get a subset of desired size

To match study requirements, we extract 100 independent samples of 100,000 individuals from the simulated population without replacement. While not fully independent, overlap is limited (~10% on average), which is deemed acceptable for the study's objectives.
## Aggregating events and exposures

Each sample is aggregated to produce death counts and central exposures :

* by age in the annuity case;

* by age and duration in the LTC case.

Aggregation is performed using the methodology described in Appendix A.

## C Derivation of marginal likelihood and LAML

This appendix presents the derivations of the marginal likelihood and its Laplace approximation (LAML) used for the automatic selection of smoothing parameters.

In the Gaussian case, where the model assumes normal conditional and prior distributions, the marginal likelihood can be computed in closed form by integrating out the latent parameters. This yields an explicit expression involving the penalty and weight matrices, and forms the basis of the outer iteration strategy.

For more general models in the exponential family, no closed-form solution is available. Instead, we apply a Laplace approximation to the marginal likelihood, based on a second-order expansion of the penalized log-likelihood around its maximum. The resulting LAML criterion is used in the performance and alternated iteration approaches for efficient and principled smoothing parameter selection.

### Marginal likelihood

Assume that $\mathbf{y} | \boldsymbol{\theta} \sim \mathcal{N}(\boldsymbol{\theta}, \sigma^2 \mathbf{W}^{-})$ and $\boldsymbol{\theta} | \lambda \sim \mathcal{N}(0, \sigma^2 \mathbf{P}_{\lambda}^{-})$. In the empirical Bayes approach, the smoothing parameter $\lambda$ is estimated by maximizing the marginal likelihood:

$$ \mathcal{L}_{\text{norm}}^{m}(\lambda) = f(\mathbf{y} | \lambda) = \int f(\mathbf{y}, \boldsymbol{\theta} | \lambda) d\boldsymbol{\theta} = \int f(\mathbf{y} | \boldsymbol{\theta}) f(\boldsymbol{\theta} | \lambda) d\boldsymbol{\theta}. $$

The conditional and prior densities involved in this integral are:

$$ f(\mathbf{y} | \boldsymbol{\theta}) = \sqrt{\frac{|\mathbf{W}|_+}{(2\pi\sigma^2)^{n_*}}} \exp \left( -\frac{1}{2\sigma^2} (\mathbf{y} - \boldsymbol{\theta})^T \mathbf{W} (\mathbf{y} - \boldsymbol{\theta}) \right) $$

$$ f(\boldsymbol{\theta} | \lambda) = \sqrt{\frac{|\mathbf{P}_{\lambda}|_+}{(2\pi\sigma^2)^{p-q}}} \exp \left( -\frac{1}{2\sigma^2} \boldsymbol{\theta}^T \mathbf{P}_{\lambda} \boldsymbol{\theta} \right) $$

where $n_*$ is the number of non-zero weights and $q$ is the order of the penalization (corresponding to the rank deficiency of $\mathbf{P}_{\lambda}$).
We then apply a second-order Taylor expansion of the joint log-density $\ln f(\mathbf{y}, \boldsymbol{\theta} | \lambda)$ around its mode $\hat{\boldsymbol{\theta}}_{\lambda}$ to approximate the integral:

$$ \ln f(\mathbf{y}, \boldsymbol{\theta} | \lambda) = \ln f(\mathbf{y}, \hat{\boldsymbol{\theta}}_{\lambda} | \lambda) + \frac{1}{2}(\boldsymbol{\theta} - \hat{\boldsymbol{\theta}}_{\lambda})^{T}(W + P_{\lambda})(\boldsymbol{\theta} - \hat{\boldsymbol{\theta}}_{\lambda}) $$

This gives the following expression for the marginal likelihood:

$$ \mathcal{L}_{\text{norm}}^{m}(\lambda) = f(\mathbf{y}, \hat{\boldsymbol{\theta}}_{\lambda} | \lambda) \int \exp \left[ -\frac{1}{2\sigma^2} (\boldsymbol{\theta} - \hat{\boldsymbol{\theta}}_{\lambda})^{T} (W + P_{\lambda}) (\boldsymbol{\theta} - \hat{\boldsymbol{\theta}}_{\lambda}) \right] d\boldsymbol{\theta} $$

which evaluates to:

$$ \mathcal{L}_{\text{norm}}^{m}(\lambda) = \sqrt{\frac{|W|_{+} |P_{\lambda}|_{+}}{(2\pi\sigma^2)^{n_{*}-q}|W + P_{\lambda}|}} \exp \left( -\frac{1}{2\sigma^2} \left[ (\mathbf{y} - \hat{\boldsymbol{\theta}}_{\lambda})^{T} W (\mathbf{y} - \hat{\boldsymbol{\theta}}_{\lambda}) + \hat{\boldsymbol{\theta}}_{\lambda}^{T} P_{\lambda} \hat{\boldsymbol{\theta}}_{\lambda} \right] \right) $$

where $\hat{\boldsymbol{\theta}}_{\lambda} = (W + P_{\lambda})^{-1}W\mathbf{y}$.

Taking the logarithm yields the marginal log-likelihood:

$$ \ell_{\text{norm}}^{m}(\lambda) = -\frac{1}{2} \left[ \frac{(\mathbf{y} - \hat{\boldsymbol{\theta}}_{\lambda})^{T} W (\mathbf{y} - \hat{\boldsymbol{\theta}}_{\lambda})}{\sigma^2} + \hat{\boldsymbol{\theta}}_{\lambda}^{T} P_{\lambda} \hat{\boldsymbol{\theta}}_{\lambda} \sigma^2 + \ln|W + P_{\lambda}| - \ln|P_{\lambda}|_{+} + C \right]. $$

where $\hat{\boldsymbol{\theta}}_{\lambda} = (W + P_{\lambda})^{-1}W\mathbf{y}$, and $C = -\ln|W|_{+} + (n_{*} - q)\ln(2\pi\sigma^2)$ is a constant independent of $\lambda$.

## Laplace approximation of the marginal likelihood

Assume that the log-likelihood $\ell(\boldsymbol{\theta})$ is combined with a Gaussian prior on the parameter vector: $\boldsymbol{\theta} \sim \mathcal{N}(0, P_{\lambda}^{-1})$. The marginal likelihood of the data ($\mathbf{d}, \mathbf{e}_{c}$) given $\lambda$ is:

$$ \mathcal{L}_{\text{ML}}^{m}(\lambda) = f(\mathbf{d}, \mathbf{e}_{c} | \lambda) = \int f(\mathbf{d}, \mathbf{e}_{c}, \boldsymbol{\theta} | \lambda) f(\boldsymbol{\theta} | \lambda) d\boldsymbol{\theta}. $$

Since no closed-form expression exists for this integral in the general exponential family case, we apply a second-order Taylor expansion of the log-posterior around its mode $\hat{\boldsymbol{\theta}}_{\lambda} = \arg\max \ell_{P}(\boldsymbol{\theta})$, where $\ell_{P}(\boldsymbol{\theta}) = \ell(\boldsymbol{\theta}) - \frac{1}{2}\boldsymbol{\theta}^{T}P_{\lambda}\boldsymbol{\theta}$ is the penalized log-likelihood.

The resulting Laplace approximation of the marginal likelihood is:

$$ \mathcal{L}_{\text{ML}}^{m}(\lambda) \approx \exp \left( \ell_{P}(\hat{\boldsymbol{\theta}}_{\lambda}) \right) \sqrt{\frac{(2\pi)^{p}}{|W_{\lambda} + P_{\lambda}|}}, $$

where $W_{\lambda} = \text{Diag}(\exp(\hat{\boldsymbol{\theta}}_{\lambda}) \odot \mathbf{e}_{c})$ is the observed Fisher information. Taking the logarithm leads to the LAML criterion:

$$ \ell_{\text{ML}}^{m}(\lambda) \approx \ell(\hat{\boldsymbol{\theta}}_{\lambda}) - \frac{1}{2} \left[ \hat{\boldsymbol{\theta}}_{\lambda}^{T} P_{\lambda} \hat{\boldsymbol{\theta}}_{\lambda} + \ln|W_{\lambda} + P_{\lambda}| - \ln|P_{\lambda}|_{+} - q \ln(2\pi) \right] \stackrel{\text{def}}{=} \ell_{\text{LAML}}^{m}(\lambda). $$
# D Algorithms

This appendix presents the computational procedures used to implement the generalized Whittaker-Henderson (WH) smoothing framework and the various automatic selection methods for the associated smoothing parameters.

## Generalized WH smoothing

Algorithm 1 implements the core iterative procedure for generalized WH smoothing, as introduced in Section 4. It details the iterative computation of the estimated log-rates $\hat{\theta}$ given fixed smoothing parameters and a chosen differencing order. The algorithm iteratively solves a penalized weighted least-squares problem until convergence is achieved, based on a predefined deviance threshold.

<table><caption>**Algorithm 1:** Iterative solution of generalized Whittaker-Henderson smoothing</caption><tr><td>inputs</td><td>: d and e<sub>c</sub></td></tr><tr><td>outputs</td><td>: ̂θ</td></tr><tr><td>parameters:</td><td>λ, q, ε<sub>d</sub> = 10<sup>-8</sup></td></tr><tr><td>begin</td><td></td></tr><tr><td>k ← 0</td><td>Construct the penalty matrix P<sub>λ</sub> based on the difference matrices of order q.</td></tr><tr><td>θ<sub>0</sub> ← ln(d*/e<sub>c</sub>), where d* = max(d, ε)</td><td></td></tr><tr><td>dev<sub>0</sub> ← ∞, cond ← true</td><td></td></tr><tr><td>while cond do</td><td></td></tr><tr><td>w<sub>k</sub> ← Diag(exp(θ<sub>k</sub>) ⊙ e<sub>c</sub>)</td><td>Form W<sub>k</sub> + P<sub>λ</sub> by adding w<sub>k</sub> to the diagonal of P<sub>λ</sub>.</td></tr><tr><td>z<sub>k</sub> ← θ<sub>k</sub> + d/w<sub>k</sub> - 1</td><td>Find the Cholesky factor R of W<sub>k</sub> + P<sub>λ</sub>.</td></tr><tr><td>Find u such that R<sup>T</sup>u = w<sub>k</sub> ⊙ z<sub>k</sub> by forward substitution.</td><td>Find u such that R<sup>T</sup>u = w<sub>k</sub> ⊙ z<sub>k</sub> by forward substitution.</td></tr><tr><td>Find θ<sub>k+1</sub> such that Rθ<sub>k+1</sub> = u by backward substitution.</td><td>Find θ<sub>k+1</sub> such that Rθ<sub>k+1</sub> = u by backward substitution.</td></tr><tr><td>dev<sub>k+1</sub> ← dev<sub>P</sub>(θ<sub>k+1</sub>), cond ← dev<sub>k+1</sub> ≤ (1 − ε<sub>d</sub>)dev<sub>k</sub></td><td>dev<sub>k+1</sub> ← dev<sub>P</sub>(θ<sub>k+1</sub>), cond ← dev<sub>k+1</sub> ≤ (1 − ε<sub>d</sub>)dev<sub>k</sub></td></tr><tr><td>k ← k + 1</td><td></td></tr><tr><td>̂θ ← θ<sub>k</sub></td><td></td></tr></table>

## Smoothing parameter selection approaches

In Section 5, we introduced three alternative strategies to automatically calibrate the smoothing parameter λ, based on marginal likelihood maximization. The following three algorithms formalize their respective procedures.
Together, these algorithms offer a modular and flexible framework for implementing WH smoothing and its data-driven calibration in practical applications.

## Outer iteration

Algorithm 2 corresponds to the outer iteration approach. In this strategy, a series of candidate values for $\lambda$ are tested sequentially. For each candidate, the generalized WH smoother of 1 is applied until convergence, and the corresponding marginal likelihood is evaluated. The process continues until no further improvement is observed. This approach separates the parameter selection and smoothing steps into nested loops.

<caption>**Algorithm 2:** Smoothing parameter selection for generalized Whittaker-Henderson smoothing - outer iteration approach.</caption>

inputs : d and e_c
outputs : $\hat{\lambda}$
parameters: $q, \epsilon_d = 10^{-8}, \epsilon_{laml} = 10^{-8}$
begin
    k ← 0
    laml_0 ← $\infty$, cond_{laml} ← true
    while cond_{laml} do
        If $k = 0$, choose an arbitrary initial value $\lambda_1$ for the smoothing parameter(s); otherwise, choose the next value $\lambda_{k+1}$ using the desired heuristic.
        k ← $k + 1$
        Use Algorithm 1 to determine the vector $\hat{\theta}_{\lambda_k}$ associated with the choice of $\lambda_k$, using a convergence threshold of $\epsilon_d$.
        Calculate the marginal likelihood $\ell_{LAML}^m(\lambda_k)$ associated with the choice of $\lambda_k$ using the intermediate quantities calculated during the estimation of $\hat{\theta}_{\lambda_k}$.
        laml_k ← $\ell_{LAML}^m(\lambda_k)$, cond_{laml} ← laml_k ≥ $(1 + \epsilon_{laml})laml_{k-1}$
    $\hat{\lambda} \leftarrow \lambda_k$

## Performance iteration

Algorithm 3 implements the performance iteration strategy. Here, the smoothing parameter is optimized at each step based on updated pseudo-response and weight vectors. The smoother and the parameter estimation are intertwined, with $\lambda$ being re-optimized after each update of the linear predictor. This often leads to faster convergence toward the maximum marginal likelihood compared to the outer iteration approach.

Let us note that in the algorithm, $\ell_{LAML}^m(\lambda_{k+1} | \theta_{k+1})$ denotes an approximate LAML, in which $\theta_{k+1}$—the maximizer of the approximate normal marginal likelihood constructed from
$\mathbf{w}_k$ and $\mathbf{z}_k$—replaces the true penalized likelihood maximizer $\hat{\theta}_{\lambda_{k+1}}$, which is not available in the performance iteration approach.

<caption>**Algorithm 3:** Parameter selection for generalized Whittaker-Henderson smoothing - performance iteration approach</caption>

**inputs** : $\mathbf{d}$ and $\mathbf{e}_c$
**outputs** : $\hat{\lambda}$
**parameters**: $q$, $\epsilon_{ml} = 10^{-8}$, $\epsilon_{laml} = 10^{-8}$
begin
    $k \leftarrow 0$
    $\theta_0 \leftarrow \ln(\mathbf{d}/\mathbf{e}_c)$
    $\text{laml}_0 \leftarrow \infty$, $\text{cond}_{\text{laml}} \leftarrow \text{true}$
    \textbf{while} $\text{cond}_{\text{laml}}$ \textbf{do}
        $\mathbf{w}_k \leftarrow \text{Diag}(\exp(\theta_k) \odot \mathbf{e}_c)$
        $\mathbf{z}_k \leftarrow \theta_k + \mathbf{d}/\mathbf{w}_k - 1$
        Find the parameter $\lambda_{k+1}$ maximizing the marginal likelihood $\ell_{\text{norm}}^m$ associated with the observation vector $\mathbf{z}_k$ and the weight vector $\mathbf{w}_k$, using the desired heuristic, using a convergence threshold of $\epsilon_{ml}$.
        Form $W_k + P_{\lambda_{k+1}}$ by adding $\mathbf{w}_k$ to the diagonal of $P_{\lambda_{k+1}}$.
        Find the Cholesky factor $R$ of $W_k + P_{\lambda_{k+1}}$.
        Find $\mathbf{u}$ such that $R^T \mathbf{u} = \mathbf{w}_k \odot \mathbf{z}_k$ by forward substitution.
        Find $\theta_{k+1}$ such that $R\theta_{k+1} = \mathbf{u}$ by backward substitution.
        $\text{laml}_{k+1} \leftarrow \ell_{\text{LAML}}^m(\lambda_{k+1} | \theta_{k+1})$, $\text{cond}_{\text{laml}} \leftarrow \text{laml}_{k+1} \ge (1 + \epsilon_{\text{laml}})\text{laml}_k$
        $k \leftarrow k + 1$
    \end{while}
    $\hat{\lambda} \leftarrow \lambda_k$;

## Alternated iteration

Finally, Algorithm 4 presents the alternated iteration approach. This hybrid method alternates between updating the smoothing parameter $\lambda$ on one hand, and an updated pseudo-response and weight vector on the other hand until convergence of an approximate marginal likelihood.

As in the performance iteration approach, $\ell_{\text{LAML}}^m(\lambda_{k+1} | \theta_{k+1})$ denotes an approximate LAML, based on the maximizer of an approximate normal marginal likelihood rather than the true penalized likelihood maximizer.

## E Illustrations for the natural parameterization of WH smoothing

This appendix provides an interpretation of Whittaker-Henderson (WH) smoothing in the framework of natural parameterization, based on the eigendecomposition of the penalty matrix.
**Algorithm 4:** Parameter selection for generalized Whittaker-Henderson smoothing -
alternated iteration approach

inputs : d and e_c

outputs : $\hat{\lambda}$

parameters: $q, \epsilon_{\text{laml}} = 10^{-8}$

begin

$k \leftarrow 0$

$\theta_0 \leftarrow \ln(d/e_c)$

$\text{laml}_0 \leftarrow \infty$, $\text{cond}_{\text{laml}} \leftarrow \text{true}$

while $\text{cond}_{\text{laml}}$ do

$\mathbf{w}_k \leftarrow \text{Diag}(\exp(\theta_k) \odot \mathbf{e}_c)$

$\mathbf{z}_k \leftarrow \theta_k + \mathbf{d}/\mathbf{w}_k - 1$

If $k = 0$, choose an arbitrary value $\lambda_1$ for the smoothing parameter(s); otherwise, choose the next value $\lambda_{k+1}$ to improve the marginal likelihood $\ell_{\text{norm}}^m$ using the desired heuristic.

Form $W_k + P_{\lambda_{k+1}}$ by adding $\mathbf{w}_k$ to the diagonal of $P_{\lambda_{k+1}}$.

Find the Cholesky factor $R$ of $W_k + P_{\lambda_{k+1}}$.

Find $\mathbf{u}$ such that $R^T \mathbf{u} = \mathbf{w}_k \odot \mathbf{z}_k$ by forward substitution.

Find $\theta_{k+1}$ such that $R\theta_{k+1} = \mathbf{u}$ by backward substitution.

$\text{laml}_{k+1} \leftarrow \ell_{\text{LAML}}^m(\lambda_{k+1} | \theta_{k+1})$, $\text{cond}_{\text{laml}} \leftarrow \text{laml}_{k+1} \ge (1 + \epsilon_{\text{laml}})\text{laml}_k$

$k \leftarrow k + 1$

$\hat{\lambda} \leftarrow \lambda_k$

This approach reveals how smoothing acts as a spectral filter that progressively attenuates
components of the signal associated with rougher variations.

In the one-dimensional case, the smoothing problem is re-expressed in a rotated basis formed
by the eigenvectors of the penalty operator, leading to a clear decomposition of the signal into
smoother and rougher components. The effect of smoothing is visualized via effective degrees
of freedom associated with each component.

The two-dimensional extension leverages Kronecker product identities to generalize the spectral
interpretation to tensor-product smoothing penalties. Illustrations are provided to highlight how
smoothing parameters affect the influence of each spectral component in both dimensions.

**One-dimensional case**

Building on the natural parameterization introduced by Demmler and Reinsch (1975), the WH
estimator can be reformulated as the solution to the following optimization problem:

$$
\hat{y} = U\hat{\beta}, \quad \hat{\beta} = \underset{\beta}{\operatorname{argmin}} \left\{ (\mathbf{y} - U\beta)^{\mathrm{T}} W (\mathbf{y} - U\beta) + \lambda \beta^{\mathrm{T}} \Sigma \beta \right\}.
$$
$\hat{\beta}$ can be interpreted as a vector of coordinates in the basis of eigenvectors of $P_{\lambda}$, yielding a spectral decomposition of the signal into components with varying degrees of smoothness, as determined by the associated eigenvalues. Figure 12 represents 8 of the eigenvectors associated with $q = 2$ for a basis of size $n = 45$. The first $q$ eigenvalues are zero as $D_{n,q}$ is of rank $n-q$.

<figure><img src="image_14.png" /><figcaption>Figure 12: A subset of the eigenvectors for the penalization matrix $D_{n,q}^T D_{n,q}$ with $n = 45$ and $q = 2$.</figcaption></figure>

By using the fact that $U^{-1} = U^T$ and by linking this expression back to the original smoothing formulation, we obtain the explicit solution:

$$ \hat{\mathbf{y}} = U(U^T W U + S_{\lambda})^{-1} U^T W \mathbf{y} \quad \text{where} \quad S_{\lambda} = \lambda \Sigma. \qquad (18) $$

In order to interpret Equation 18, consider the special case where all weights are equal to 1 and therefore:

$$ \hat{\mathbf{y}} = U(U^T U + \lambda \Sigma)^{-1} U^T \mathbf{y} = U(I_n + \lambda \Sigma)^{-1} U^T \mathbf{y}. $$

The transformation from $\mathbf{y}$ to $\hat{\mathbf{y}}$ can then be seen as a 3-step process, reading the equation from right to left:

1. Decomposition of the signal $\mathbf{y}$ in the basis of eigenvectors through the left multiplication by $U^T$.
2. Attenuation of the signal components based on the eigenvalues associated with these components. If we denote $s = \text{diag}(\Sigma)$, then $(I_n + \lambda\Sigma)^{-1} = \text{Diag}[1/(1 + \lambda s)]$. After the left multiplication by $(I_n + \lambda\Sigma)^{-1}$, each component is hence divided by a factor $1 + \lambda s \ge 1$. This coefficient increases linearly with $\lambda$, but the rate of increase varies with the magnitude of the corresponding eigenvalue.

3. Recomposition of the attenuated signal in the canonical basis through the left multiplication by $U$.

When weights are not uniform, the structure becomes more complex since $U^T W U$ is no longer a diagonal matrix. However, it is still possible to interpret the effect of smoothing from the diagonal of the matrix $F = (U^T W U + S_\lambda)^{-1} U^T W U$. Indeed:

$$U^T \hat{\mathbf{y}} = U^T U \hat{\mathbf{\theta}} = \hat{\mathbf{\theta}} = (U^T W U + S_\lambda)^{-1} U^T W \mathbf{y} = (U^T W U + S_\lambda)^{-1} U^T W U U^T \mathbf{y} = F U^T \mathbf{y}.$$

Since the vectors $U^T \mathbf{y}$ and $\hat{\mathbf{\theta}} = U^T \hat{\mathbf{y}}$ represent the coordinates of $\mathbf{y}$ and $\hat{\mathbf{y}}$ respectively in the basis of eigenvectors of $D_{n,q}^T D_{n,q}$, $F$ thus acts as a transformation matrix on the spectral coordinates, analogous to the role played by the hat matrix $H = U (U^T W U + S_\lambda)^{-1} U^T W$ for the observations. The diagonal values of $F$ may be interpreted as the effective degrees of freedom associated with each eigenvector after smoothing. It can be verified that:

$$\text{tr}(F) = \text{tr}[(U^T W U + S_\lambda)^{-1} U^T W U] = \text{tr}[U (U^T W U + S_\lambda)^{-1} U^T W] = \text{tr}(H)$$

which means that the sum of the effective degrees of freedom remains the same whether it is counted per observation or per parameter.

Figure 13 represents the effective degrees of freedom per parameter in the previous illustration of smoothing. The first $q$ eigenvectors are never penalized, so their effective degrees of freedom are always equal to 1, regardless of the smoothing parameter used. The other eigenvectors have strictly decreasing effective degrees of freedom with $\lambda$. These degrees of freedom are mostly decreasing with increasing eigenvalues of $D_{n,q}^T D_{n,q}$, although in the presence of non-unit weights and for small values of $\lambda$, this is not always the case.

## Two-dimensional case

In the two-dimensional case, we have $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}$. Similar to the one-dimensional case, we can perform the eigendecomposition of the matrices $D_{n_x,q_x}^T D_{n_x,q_x}$ and $D_{n_z,q_z}^T D_{n_z,q_z}$, yielding $D_{n_x,q_x}^T D_{n_x,q_x} = U_x \Sigma_x U_x^T$ and $D_{n_z,q_z}^T D_{n_z,q_z} = U_z \Sigma_z U_z^T$. Define $U = U_z \otimes U_x$ and perform the reparametrization $\beta = U^T \mathbf{\theta} \Leftrightarrow \mathbf{\theta} = U \beta$. By leveraging the properties of the Kronecker product, we can rewrite the smoothness criterion in a simplified form:

$$\mathbf{\theta}^T P_\lambda \mathbf{\theta} = (U \beta)^T P_\lambda (U \beta) = \beta^T (\lambda_x I_{n_z} \otimes \Sigma_x + \lambda_z \Sigma_z \otimes I_{n_x}) \beta.$$

This leads to an alternative formulation of the optimization problem:
<figure><img src="image_15.png" /><figcaption>Figure 13: Effective degrees of freedom per eigenvector after applying 1D WH smoothing, for different combinations of smoothing parameter.</figcaption></figure>

$$ \hat{y} = U\hat{\beta}, \quad \hat{\beta} = \underset{\beta}{\operatorname{argmin}} \left\{ (\mathbf{y} - U\beta)^{\mathrm{T}} W (\mathbf{y} - U\beta) + \lambda \beta^{\mathrm{T}} (\lambda_{x} I_{n_{x}} \otimes \Sigma_{x} + \lambda_{z} \Sigma_{z} \otimes I_{n_{z}}) \beta \right\}. $$

The solution to the smoothing problem, as in the one-dimensional case, is given by:

$$ \hat{y} = U(U^{\mathrm{T}}WU + S_{\lambda})^{-1}U^{\mathrm{T}}Wy \quad \text{where} \quad S_{\lambda} = \lambda_{x}I_{n_{x}} \otimes \Sigma_{x} + \lambda_{z}\Sigma_{z} \otimes I_{n_{z}}. $$

Figure 14 represents the residual degrees of freedom associated with each parameter after applying the smoothing, in the two-dimensional case, for different combinations of the smoothing parameters. Similar to the one-dimensional case, these degrees of freedom decrease as the smoothing parameters increase and are particularly small for higher eigenvalues. The eigenvectors are sorted in ascending order of eigenvalues for each one-dimensional penalty matrix $D_{n_x,q_x}^T D_{n_x,q_x}$ and $D_{n_z,q_z}^T D_{n_z,q_z}$.

## F Derivation of constrained extrapolation in the 2D case

This appendix provides a closed-form expression for the constrained extrapolated estimator, which extends the smoothed surface beyond the observed domain while preserving the values estimated during the initial fit. It also derives the associated variance-covariance matrix, accounting for both propagated uncertainty and additional variability in the extrapolated region.

To obtain an estimator $\hat{y}_+$ that minimizes the extended optimization problem under the constraint of preserving the initial coefficients, i.e. $C\hat{y}_+^* = \hat{y}$, we follow the approach proposed by Carballo, Durban, and Lee (2021) and introduce a Lagrange multiplier $\omega$. The associated constrained extended optimization problem is now written as:
<figure><img src="image_16.png" /><figcaption>Figure 14: Residual degrees of freedom per eigenvector after applying 2D WH smoothing, for different combinations of smoothing parameters.</figcaption></figure>

$$ (\hat{\mathbf{y}}_+^*, \hat{\omega}) = \operatorname*{argmin}_{\mathbf{\theta}_+^*, \omega} \left\{ (\mathbf{y}_+ - \mathbf{\theta}_+^*)^T W_+(\mathbf{y}_+ - \mathbf{\theta}_+^*) + \mathbf{\theta}_+^{*T} P_+ \mathbf{\theta}_+^* + 2\omega^T (C \mathbf{\theta}_+^* - \hat{\mathbf{y}}) \right\}. \quad (19) $$

Taking the partial derivatives of Equation 19 with respect to $\mathbf{\theta}_+^*$ and $\omega$ gives:

$$ \frac{\partial}{\partial \mathbf{\theta}_+^*} \left\{ (\mathbf{y}_+ - \mathbf{\theta}_+^*)^T W_+(\mathbf{y}_+ - \mathbf{\theta}_+^*) + \mathbf{\theta}_+^{*T} P_+ \mathbf{\theta}_+^* + 2\omega^T (C \mathbf{\theta}_+^* - \hat{\mathbf{y}}) \right\} = -2W_+(\mathbf{y}_+ - \mathbf{\theta}_+^*) + 2P_+ \mathbf{\theta}_+^* + 2\omega^T C $$

$$ \frac{\partial}{\partial \omega} \left\{ (\mathbf{y}_+ - \mathbf{\theta}_+^*)^T W_+(\mathbf{y}_+ - \mathbf{\theta}_+^*) + \mathbf{\theta}_+^{*T} P_+ \mathbf{\theta}_+^* + 2\omega^T (C \mathbf{\theta}_+^* - \hat{\mathbf{y}}) \right\} = 2(C \mathbf{\theta}_+^* - \hat{\mathbf{y}}) $$

Setting these derivatives to zero yields the linear system:

$$ \begin{bmatrix} W_+ + P_+ & C^T \\ C & 0 \end{bmatrix} \begin{bmatrix} \hat{\mathbf{y}}_+^* \\ \hat{\omega} \end{bmatrix} = \begin{bmatrix} W_+ \mathbf{y}_+ \\ \hat{\mathbf{y}} \end{bmatrix} $$
The solution for $\hat{\mathbf{y}}_+^*$ can be derived using formulas for the inversion of a symmetric matrix partitioned with 2 × 2 blocks, setting $V_+ = (W_+ + P_+)^{-1}$:

$$ \hat{\mathbf{y}}_+^* = V_+ \{ I - C^T [CV_+ C^T]^{-1} CV_+ \} W_+ \mathbf{y}_+ + V_+ C^T [CV_+ C^T]^{-1} \hat{\mathbf{y}}. $$

Since $W_+ = C^TWC$, the first term is actually zero, and this expression simplifies to:

$$ \hat{\mathbf{y}}_+^* = V_+ C^T [CV_+ C^T]^{-1} \hat{\mathbf{y}} = V_+ C^T (V_+^{11})^{-1} \hat{\mathbf{y}} = Q^T \left[ -(P_+^{22})^{-1} P_+^{21} \right] \hat{\mathbf{y}} $$

which is a linear transformation of $\hat{\mathbf{y}}$.

Defining $A_+^* = Q^T \left[ \begin{matrix} I \\ -(P_+^{22})^{-1} P_+^{21} \end{matrix} \right]$, a variance-covariance of $\mathbf{y}_+^* | \theta_+$ based on this expression is given by:

$$ A_+^* V A_+^{*T} = Q^T \left[ \begin{matrix} V & -V P_+^{12} (P_+^{22})^{-1} \\ -(P_+^{22})^{-1} P_+^{21} V & (P_+^{22})^{-1} P_+^{21} V P_+^{12} (P_+^{22})^{-1} \end{matrix} \right] Q. \quad (20) $$

Equation 20 is very similar to the equivalent formulation in the one-dimensional case with however two main differences. First, every occurrence of $V_+^{11}$ is replaced by $V$. This is consistent with the constraint that the coefficients over the initial domain are held fixed at their estimated values. Second, as the solution to the constrained extended optimization problem of Equation 19 was expressed as a linear transformation of $\hat{\mathbf{y}}$, Equation 20 is missing the innovation error term $(P_+^{22})^{-1}$ associated with the prior on the extrapolated coefficients. Not including this term would be tantamount to considering that $\theta_+$ has some degree of variability in the region of the initial data but is perfectly smooth beyond this range. Adding this innovation error, we obtain the following variance-covariance matrix for the constrained optimization problem:

$$ V_+^* = Q^T \left[ \begin{matrix} V & -V P_+^{12} (P_+^{22})^{-1} \\ -(P_+^{22})^{-1} P_+^{21} V & (P_+^{22})^{-1} P_+^{21} V P_+^{12} (P_+^{22})^{-1} + (P_+^{22})^{-1} \end{matrix} \right] Q. $$

which still verifies

$$ V_+^* W_+ \mathbf{y}_+ = V_+^* C W \mathbf{y} = Q^T \left[ \begin{matrix} I \\ -(P_+^{22})^{-1} P_+^{21} \end{matrix} \right] V W \mathbf{y} = \hat{\mathbf{y}}_+^*. $$

## G Choosing the order of the difference matrices

This appendix investigates how the choice of the difference order $q$ in WH smoothing affects model performance. It provides empirical guidance for selecting $q$ based on AIC values computed across simulation replicates, and confirms that second-order differences offer a good balance between fit quality and extrapolation stability.
In Whittaker-Henderson (WH) smoothing, the order of the difference matrix used in the penalization term, typically denoted $q$ (or $(q_x, q_z)$ in the two-dimensional case), governs the smoothness prior imposed on the underlying signal. While Whittaker originally proposed using third-order differences, second-order differences have since become the standard choice. This aligns with the common statistical definition of smoothness as the integrated squared second derivative, a justification also used for cubic spline smoothing (Reinsch 1967) and P-splines (Eilers and Marx 1996).

From a Bayesian perspective, penalizing second-order differences is equivalent to assuming that the log-transformed target function follows a locally affine (i.e., linear) trajectory. For mortality data, this implies a prior belief in exponential growth of the mortality rate—a widely accepted assumption consistent with the Gompertz model. For other biometric risks, such as disability or long-term care, this prior is also appropriate for the age dimension if included.

Penalizations of order $q > 3$ are harder to justify in actuarial practice. They may produce undesirable artefacts in the extrapolated regions, where the extrapolation from the fit tends to follow a polynomial of degree $q - 1$.

It is possible to treat the choice of $q$ as a model selection problem. REML is not suitable in this context, as it assumes an equal fixed effect structure, however standard criteria such as AIC and GCV can be used instead. Relying on AIC and the 100 replicate datasets described in Section 6, we found that:

* For the annuity datasets, Figure 15 shows that second-order differences ($q = 2$) provides the best results.

* For the duration dimension in LTC datasets, Figure 16 shows that, regarding the age dimension, second-order differences ($q = 2$) is still the best option while for the duration dimension, AIC improves significantly when moving from first- to second-order penalization, and then slightly improves further for higher orders. This suggests that higher-order prior may be relevant for the duration dimension in this particular dataset.

Still, the gains from increasing the order beyond 2 are limited and may not outweigh the risks of erratic extrapolation. Therefore, second-order differences remain a pragmatic and robust default.
<figure><img src="image_17.png" /><figcaption>Figure 15: AIC values for different penalization orders in WH smoothing, based on 100 replicates of the annuity portfolio (100,000 individuals).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 16: AIC values for combinations of difference matrix orders along age and duration in WH smoothing, based on 100 replicates of the LTC portfolio (100,000 individuals).</figcaption></figure>