scaletikzpicturetowidth[1]

# Catastrophe risk in a stochastic multi-population mortality model

Jens Robben  Faculty of Economics and Business, KU Leuven, Belgium.  LRisk,
Leuven Research Center on Insurance and Financial Risk Analysis, KU Leuven,
Belgium.  Corresponding author<EMAIL_ADDRESS>Katrien Antonio
Faculty of Economics and Business, KU Leuven, Belgium.  Faculty of Economics
and Business, University of Amsterdam, The Netherlands.  LRisk, Leuven
Research Center on Insurance and Financial Risk Analysis, KU Leuven, Belgium.
LStat, Leuven Statistics Research Center, KU Leuven, Belgium.  RCLR, Research
Centre for Longevity Risk, University of Amsterdam, The Netherlands.

###### Abstract

This paper incorporates mortality shocks in the scenarios for future mortality
rates produced by a stochastic multi-population mortality model. Hereto, the
proposed model combines a decreasing stochastic mortality trend with a
mechanism that switches between regimes of low and high volatility. During the
high volatility regimes, mortality shocks occur that last from one to several
years and temporarily impact the mortality rates before returning to the
overall mortality trend. Furthermore, we account for the age-specific impact
of these mortality shocks on mortality rates. Actuaries and risk managers can
tailor this scenario generator to their specific needs, risk management
objectives, or supervisory requirements. The generated scenarios allow
(re)insurers, policymakers, or actuaries to evaluate the effects of different
catastrophe risk scenarios on, e.g., the calculation of solvency capital
requirements.

##### JEL classification:

G22

##### Keywords:

mortality shocks; regime switch; multi-population mortality model; mortality
improvement model; stochastic mortality model; catastrophe risk; scenario
analysis

##### Data and code availability statement:

The data sets used in this paper are publicly available from the Human
Mortality Database (https://www.mortality.org/) and Eurostat
(https://ec.europa.eu/eurostat). The code used for the implementation and
analysis of the case study presented in this paper is available on the GitHub
repository https://github.com/RobbenJ/catastrophe-risk. The code is written in
R and can be accessed and downloaded for further reference and replication of
the obtained results.

##### Funding statement:

The authors acknowledge funding from the FWO and F.R.S.-FNRS under the
Excellence of Science (EOS) programme, project ASTeRISK (40007517). The FWO
WOG network W001021N is acknowledged as well. Katrien Antonio acknowledges
support from the Chair of Excellence on Digital Insurance And Long‑term risk
(DIALog) by CNP Assurances and from the RCLR, the Research Centre for
Longevity Risk at the University of Amsterdam.

##### Conflict of interest disclosure:

The authors declare no conflict of interest.

##  1 Introduction

The overall evolution of mortality rates over time is referred to as the
mortality trend. Since the beginning of the 19th century, mortality rates have
declined in Europe, but occasional and temporary mortality jumps or shocks
have disrupted this trend (Lee, 2003). Figure 1 shows the age-averaged death
rates by aggregating total deaths and population exposures over age ranges
20-59 and 60-85 for a set of West and North European countries. The figure
highlights the main mortality shocks that happened throughout the historical
period 1850-2021. The timeline in Figure 2 reveals the nature of these shocks:
some of the observed historical shocks are related to conflicts like World War
I and World War II, while others are related to pandemics and epidemics such
as a cholera outbreak, the Spanish flu, or, more recently, the COVID-19
pandemic. As Figure 1 indicates, these shocks can be seen as sudden, upward
spikes in the mortality rates, potentially spanning multiple years, followed
by a subsequent return to the overall declining mortality trend.

Accurately modeling future mortality rates is crucial for the valuation of
life-contingent risks. To safeguard the solvability within the European
(re)insurance industry, the solvency capital requirements (SCRs) in the
Solvency II directive require the calculation of a separate SCR for mortality
risk, longevity risk, and catastrophe risk. Mortality risk refers to the risk
of deterioration in mortality rates, longevity risk is the risk of people
living longer than expected, and catastrophe risk captures sudden extreme
events that lead to a sharp increase in mortality rates, such as a nuclear
explosion or pandemic. The main objective of this paper is the development of
a toolbox that allows (re)insurers to generate scenarios for future mortality
rates, equipped with shocks that can impact the mortality trend over several
years while taking into account the age-specific impact of such a shock on
mortality rates. This enables the valuation of risks in the presence of
catastrophic shocks. Scenarios involving abrupt improvements, such as a sudden
cure for diseases like cancer (Gielen and de Waegenaere, 2014), are not
considered in this paper.

Figure 1: The age-averaged death rates, defined as
$1/|\mathcal{X}|\sum_{x\in\mathcal{X}}d_{x,t}^{A}/E_{x,t}^{A}$, with
$\mathcal{X}$ the age range under consideration, $d_{x,t}^{A}$ the total
deaths, and $E_{x,t}^{A}$ the population exposure, aggregated over a set of
West and North European countries, i.e. Austria, Belgium, Denmark, Germany,
Finland, France, Ireland, Iceland, Luxembourg, the Netherlands, Norway,
Sweden, Switzerland, and the United Kingdom. We average the death rates over
the age range 20-59 (left panel) and 60-85 (right panel). The color scheme
used for highlighting the main mortality shocks is also used in Figure 2 to
visualize the occurrence of the shocks over time. 1854-1855: Crimean war
cholera outbreak 1859: Cholera outbreak 1866: Cholera outbreak 1870-1871:
Franco-German War Smallpox epidemic 1889-1892: Asian flu 1914-1919: World War
I Spanish flu 1940-1945: World War II 2020-2021: COVID-19 pandemic Figure 2:
Historical mortality shocks in Europe from 1850-2021.

Many insurers and reinsurers worldwide express concerns regarding the
potential threats posed by catastrophic mortality events. Toole (2007) finds
that under an extreme pandemic scenario, the US life insurance industry would
face additional claims expenses up to $25\%$ of its risk-based capital, which
could cause insolvency in the (re)insurance industry. Therefore, (re)insurers
might consider, in an attempt to mitigate their exposure, investing in so-
called mortality-linked securities, e.g., the Swiss Re CAT bond (2003).
Pricing such mortality-linked securities requires insight into the modeling
and projection of mortality rates. The inclusion of mortality shocks into
mortality projection models is of significant importance. Chen and Cox (2009)
provide a case study on pricing mortality-linked securities with their
proposed mortality model, including a shock component. Zhou et al. (2013) use
an economic-pricing framework to investigate how mortality shocks can affect
the supply and demand of mortality-linked securities. Moreover, Liu and Li
(2015) provide an application for pricing hypothetical catastrophic mortality
bonds.

The proposed mortality modeling framework consists of two building blocks. The
first building block captures the overall mortality trend using a stochastic
multi-population model. The second building block accounts for catastrophe
risk via mortality shocks. We utilize the projection model equipped with
mortality shocks to generate scenarios for future mortality rates that are
exposed to the occurrence of shocks. Our model introduces shocks that can last
from one to several years, temporarily impacting the mortality rates before
returning to the overall mortality trend. Furthermore, we account for the age-
specific impact of a mortality shock on mortality rates. As such, we enable
actuaries and risk managers to tailor the occurrence and impact of these
shocks to their specific needs, risk appetite, or supervisory requirements.
This flexibility builds on the calibration of the model, where users can make
informed decisions about the inclusion or exclusion of historical mortality
shocks. For instance, users may choose to disregard the mortality shocks
caused by both world wars or assign more weight to recent observations. Such
modifications directly influence the frequency and severity of future
mortality shocks in the model projection step, empowering users to tailor the
framework to their risk management objectives.

The use of stochastic single-population mortality models originated from the
seminal work of Lee and Carter (1992), with extensions proposed by (among
others) Renshaw and Haberman (2003), Currie (2006), Renshaw and Haberman
(2006), Cairns et al. (2009), Plat (2009), Haberman and Renshaw (2011),
Levantesi and Pizzorusso (2019) and Dowd et al. (2020). In addition to single
population models, incorporating mortality data from multiple populations with
similar characteristics has been proposed as an effective way to enhance the
stability of mortality projections. While it is generally expected that
variations in mortality rates between such populations will stabilize over
time, unforeseen factors or events can still lead to temporary divergences.
The importance of multi-population models has been discussed in Antonio et al.
(2017) among others. A widely used multi-population model is the one proposed
by Li and Lee (2005). The Dutch Actuarial Association and the Institute of
Actuaries in Belgium have both adopted a stochastic multi-population mortality
model of type Li-Lee, as reported in Royal Dutch Actuarial Association (2022)
and Antonio et al. (2020) respectively. Haberman et al. (2014) discuss a
variety of other multi-population mortality models as extensions of commonly
used single-population models. Stochastic mortality models usually incorporate
a series of age and period-specific parameters to model the logarithm of the
age and period-specific force of mortality. This paper will refer to these
parameters as age and period effects. Subsequently, time series models are
used to extrapolate the calibrated period effects, which enables the
projection of future mortality rates.

We use a mortality improvement model to examine the changes in mortality rates
over time. This is in contrast to a classic mortality model that directly
specifies the mortality rate or force of mortality. We motivate this set-up by
previous work from Mitchell et al. (2013) and Hunt and Villegas (2021).
Mitchell et al. (2013) conduct an extensive empirical study comparing their
proposed improvement model with traditional mortality models for the logarithm
of the central death rate such as the Lee-Carter model (Lee and Carter, 1992),
(variants of) the model proposed by Renshaw and Haberman (2006) and the
mortality model of Plat (2009). Their study covers a comprehensive dataset
spanning 100 years across 11 countries. The results of their investigation
demonstrate the superiority of the improvement model over the aforementioned
models using two in-sample goodness-of-fit measures, i.e. the root sum of
squared errors and the unexplained variance.111The unexplained variance at age
(group) $x$ represents the proportion of variability in the observed mortality
rates at age $x$ that is not accounted for by the calibrated model. Hunt and
Villegas (2021) argue that practitioners, such as life insurers and pension
funds, are primarily interested in mortality improvement rates, especially
when assessing longevity risk. A better-than-average mortality improvement
rate that persists for several years indicates that life expectancy may exceed
initial projections. As a result, reserves may have to increase to safeguard
the solvability of the company. In the UK, the Continuous Mortality
Investigation (CMI) introduced the concept of mortality improvement rates in
2002 (Continuous Mortality Investigation Bureau, 2002) and used it in their
subsequent mortality projections. Furthermore, in the US, the Society of
Actuaries’ Retirement Plans Experience Committee has published since the year
2014 annual mortality improvement scales for projecting mortality rates.222The
Mortality Improvement Scale MP-2014 can be consulted on
https://www.soa.org/resources/experience-studies/2014/research-2014-mp/, where
the mortality improvement model and corresponding assumptions are discussed.
This shows the interest in as well as practical relevance of mortality
improvement models.

A limited body of literature explores the incorporation of mortality shocks
into stochastic mortality models. A first stream of research primarily puts
focus on detecting outliers in the time series for the calibrated period
effects. These outliers, or mortality shocks, are considered non-repetitive,
exogenous intervention effects that are not accounted for in future mortality
projections. For instance, Lee and Carter (1992) utilize a random walk with
drift in conjunction with an intervention term to capture the 1918 influenza
epidemic. Li and Chan (2005) extend this work and unravel the underlying
mortality trend by identifying different types of outliers in the residuals of
an ARIMA($p$,$d$,$q$) model fitted to the calibrated period effects. Their
outlier detection strategy identifies different types of outliers. A second
stream of literature argues that mortality shocks, including pandemic shocks,
should not be seen as exogenous, one-time events unlikely to occur again in
the future. Rather, they should be considered as recurrent events and
explicitly incorporated into mortality projection models. For this purpose,
Cox et al. (2006) and Chen and Cox (2009) augment the Lee-Carter framework
with jump effects. The magnitude of these jump effects is modeled with
independent and normally distributed random variables, while their occurrence
is determined by a Bernoulli-distributed random variable. In Cox et al.
(2006), the jump effect is considered to be permanent as the magnitude of a
particular jump impacts all future mortality rates. However, Chen and Cox
(2009) argue that most mortality jumps or shocks are transitory in nature and
implement the jump effect accordingly, affecting the mortality rates only
during a single year. Zhou et al. (2013) further extend the work of Chen and
Cox (2009) to a two-population framework. Since these jump effects are
included in the time series model for the period effect(s), Cox et al. (2006),
Chen and Cox (2009), and Zhou et al. (2013) implicitly assume in their work
that the age effect of a mortality shock is identical to the age effect
observed in the general mortality trend of the Lee-Carter model. Liu and Li
(2015) point out that this is a too restrictive assumption that does not align
with empirical observations. Their analysis reveals that mortality shocks in
the United States and England and Wales exhibit varying age patterns,
deviating from the age effect calibrated with a Lee-Carter model. Furthermore,
distinct mortality shocks also tend to impact different age groups, e.g., a
war mainly affects younger ages, whereas the recent COVID-19 pandemic
particularly affected older ages. Hence, Liu and Li (2015) extend the Lee-
Carter modeling framework with a jump effect and dedicated age-specific
parameters. Goes et al. (2023) introduce a Bayesian generalization of the
model proposed by Liu and Li (2015) and integrate a vanishing jump effect on
the mortality rates. Given the scarcity of severe mortality improvements and
shocks, a third stream of literature builds on the principles of extreme value
theory. Chen and Cummins (2010) and Gungah and Narsoo (2021) use a generalized
Pareto distribution to characterize the first-order differences of the Lee-
Carter’s calibrated period effect above a specific threshold, while a random
walk with drift is used below the threshold. Hence, these models effectively
capture both extreme longevity improvements as well as mortality shocks.
Furthermore, by taking random draws from the fitted generalized Pareto
distribution, the model accounts for extreme mortality events in the future.
The research avenues discussed so far fail in situations involving a prolonged
calibration period during which mortality rates may be impacted by shocks as
well as shifts in the mortality trend. In this regard, a fourth stream of
research focuses on the use of regime-switching models for the calibrated
period effects. Milidonis et al. (2011) employ a regime-switching model to
capture regimes in the mortality trend for the US population (1901-1999).
Hereto, the calibrated period effect of a Lee-Carter model is governed by a
regime-switching random walk model, where the volatility of the random walk’s
error term varies across regimes. Similarly, Hainaut (2012) proposes a regime-
switching random walk model for the calibrated period effects of a multi-
factor Lee-Carter model for the French population (1946-2007). Here, the drift
of the random walk and the volatility of the error term vary across regimes.

Our paper contributes to this literature in four ways. First, we directly
incorporate the regime switch in the mortality model specifications rather
than in the calibrated period effects, as is done by Milidonis et al. (2011)
and Hainaut (2012). To the best of our knowledge, this has not been tackled
yet in the current literature. In this way, we aim to construct a mortality
model that switches between a low volatility regime in which no mortality
shocks occur and a high volatility regime in which mortality shocks take
place. In contrast to Milidonis et al. (2011) and Hainaut (2012), we model at
the same time a dedicated set of age-specific effects that describe the impact
of the mortality shocks on the mortality rates. The latter is in line with the
extension of Liu and Li (2015) to the work of Cox et al. (2006) and Chen and
Cox (2009). However, we now apply this idea to regime-switching models rather
than jump models to easily switch between periods of high and low volatility.
We further add geometrically decaying weights in the calibration process of
our mortality model to put more weight on the most recent mortality
improvement rates. This enables us to capture and extrapolate the dynamic
nature of mortality trends more effectively. Second, we calibrate our proposed
model in a multi-population setting with a substantially long calibration
period, starting from the year 1850 onwards. We use a multi-factor stochastic
multi-population mortality improvement model to capture the dynamics over this
lengthy calibration period. Third, we hereby derive an invariant set of
parameter transformations in a two-factor mortality improvement model such
that the transformed parameters satisfy the imposed identifiability
constraints, inspired by Hunt and Blake (2020), and lead to the same model
fit. Fourth, in a case study, we shed light on the calculation of the SCR in
the Solvency II directive for longevity, mortality, and catastrophe risk with
both the standard model as well as our proposed mortality model equipped with
mortality shocks.

This paper is organized as follows: Section 2 introduces the notation. Section
3 presents the proposed mortality improvement model, which is composed of a
baseline mortality improvement model and a regime-switching model with age-
specific impact. We explain the methodology for calibrating the baseline model
in Section 4.1, while Section 4.2 covers the calibration of the regime-
switching model. Section 5 explains how projections can be made with the
calibrated mortality improvement model. In the case study discussed in Section
6, we model and project Dutch mortality rates based on the proposed multi-
population mortality model equipped with mortality shocks. Furthermore, we
compare the SCR calculated using our proposed mortality improvement model with
the SCR obtained with the standard model.

##  2 Notations

Let $d_{x,t}^{(c)}$ denote the observed number of deaths in a country $c$ at
age $x$ in year $t$.333The data on deaths can be categorized as male, female,
or total (unisex) deaths. Although our case study in Section 6 concentrates
solely on male data, it is important to note that the methodology is also
applicable to both female and unisex data. Therefore, we will omit the
dependence on gender in all notations used throughout the paper. We consider a
set of possible ages $x\in\mathcal{X}$ and a range of years $t\in\mathcal{T}$.
The observed exposure in country $c$, denoted by $E_{x,t}^{(c)}$, reflects the
total amount of person-years lived by individuals aged $[x,x+1)$ through year
$[t,t+1)$. We denote the country-specific force of mortality by
$\mu_{x,t}^{(c)}$ and the central death rate by $m_{x,t}^{(c)}$. The empirical
estimate of the central death rate, also called the crude central death rate,
at age $x$ in year $t$ is defined as:

$\hat{m}_{x,t}^{(c)}=\dfrac{d_{x,t}^{(c)}}{E_{x,t}^{(c)}}.$

Furthermore, we define $q_{x,t}^{(c)}$ as the mortality rate at exact age $x$
in year $t$, i.e. the probability that an individual from country $c$, born on
January 1 of year $t-x$, who is alive at January 1 of year $t$ will die within
the next year. In this paper, we assume that the number of deaths random
variable $D_{x,t}^{(c)}$ follows a Poisson distribution (Brouhns et al.,
2002):

$\displaystyle
D_{x,t}^{(c)}\sim\text{Pois}\left(E_{x,t}^{(c)}\cdot\mu_{x,t}^{(c)}\right).$
(2.1)

The force of mortality is assumed to be constant over each age $x$ and year
$t$, i.e. $\mu_{x+s,t+s}^{(c)}=\mu_{x,s}^{(c)}$ for $s\in[0,1)$. Under these
assumptions, the maximum likelihood estimate of the force of mortality
coincides with the crude central death rate, i.e.
$\hat{\mu}_{x,t}^{(c)}=\hat{m}_{x,t}^{(c)}$ (Pitacco et al., 2009).
Additionally, the mortality rates can then be approximated by:

$\displaystyle\hat{q}_{x,t}^{(c)}=1-\exp\left(-\hat{\mu}_{x,t}^{(c)}\right).$
(2.2)

Single population stochastic mortality models typically structure the country-
specific force of mortality $\mu_{x,t}^{(c)}$ using a combination of age-
specific effects $\beta_{x}^{(j,c)}$, period-specific effects
$\kappa_{t}^{(j,c)}$ and, occasionally, cohort-specific effects
$\gamma_{t-x}^{(j,c)}$, for $j\in\\{1,...,l\\}$, as:

$\displaystyle\log\mu_{x,t}^{(c)}=\beta_{x}^{(1,c)}\kappa_{t}^{(1,c)}\gamma_{t-x}^{(1,c)}+\ldots+\beta_{x}^{(l,c)}\kappa_{t}^{(l,c)}\gamma_{t-x}^{(l,c)},$
(2.3)

where, potentially, some age, period, or cohort effects are fixed to zero or
one. A common approach for estimating the involved parameters is to maximize
the Poisson log-likelihood derived from Equation (2.1). A stochastic mortality
model then typically specifies and calibrates time series models to project
the estimated period and cohort-specific effects.

Multi-population stochastic mortality models are widely used to examine
mortality patterns across populations that share similar characteristics,
offering enhanced stability in mortality projections. These models structure
the country-specific force of mortality by combining a set of country-specific
age, period, and cohort effects, as in Equation (2.3), with a set of common
age, period, and cohort effects. The latter structure a common, multi-
population mortality trend that is shared by all the populations under
consideration, whereas the country-specific parameters capture deviations from
this common trend for the particular country of interest.

Recently, more attention has been given to stochastic mortality improvement
models. Such improvement models structure the (force of) mortality improvement
rates, defined as:

$\displaystyle\log\left(\frac{\mu_{x,t}^{(c)}}{\mu_{x,t-1}^{(c)}}\right)=\log\mu_{x,t}^{(c)}-\log\mu_{x,t-1}^{(c)}.$
(2.4)

A positive value for the mortality improvement rate at year $t$ indicates a
mortality deterioration relative to the preceding year. Conversely, a negative
value indicates an improvement in mortality compared to year $t-1$. Mortality
improvement models therefore aim to capture the variability in the rate of
change in mortality over time, as opposed to the level of mortality at a
specific time point.

##  3 A stochastic multi-population mortality improvement model with a shock
regime

The proposed mortality model comprises a baseline mortality improvement model
and a regime-switching model, as pictured in Figure 3. In the first stage, we
design and estimate the baseline mortality improvement model to capture the
country-specific overall mortality decline, which is illustrated on the left-
hand side of the figure. The residuals obtained from this estimated baseline
improvement model are vulnerable to periods of high and low volatility, as
evidenced on the right-hand side of Figure 3. The high volatility periods
arise due to country-specific mortality shocks caused by war or pandemic
events. Therefore, in the second stage, we use a regime-switching model to
capture such shocks present in the residuals. Hereto, an underlying Markov
chain alternates between a low volatility state (LVS) and a high volatility
state (HVS). The LVS mode tracks the country-specific mortality decline, while
the HVS mode accounts for the shocks.

$LVS$ $HVS$ $1-p_{12}$$p_{12}$$1-p_{21}$$p_{21}$$t$

Figure 3: Visual representation of the proposed mortality improvement model,
consisting of a baseline mortality improvement model (left) and a regime-
switching model (right). The baseline improvement model captures the country-
specific mortality decline. The regime-switching model captures the country-
specific mortality shocks and is constructed on the residuals of the estimated
baseline model. In this way, the mortality improvement rates alternate between
a low volatility state (LVS), which tracks the country-specific mortality
decline, and a high volatility state (HVS), which captures the occurrence of
mortality shocks.

###  3.1 The baseline mortality improvement model

##### Model specification.

We introduce our mortality improvement model as an extension of a multi-factor
multi-population model (Renshaw and Haberman, 2003):

$\displaystyle\log\mu_{x,t}^{(c)}-\log\mu_{x,t-1}^{(c)}=A_{x}+\displaystyle\sum_{i=1}^{m}B_{x}^{(i)}K_{t}^{(i)}+\displaystyle\sum_{j=1}^{l}\beta_{x}^{(j,c)}\kappa_{t}^{(j,c)}.$
(3.1)

We use a set of $m$ age/period effects to describe a common, multi-population
mortality improvement trend. Subsequently, a set of $l$ age/period effects
captures country-specific deviations from this common improvement trend.
Depending on the area of interest, other types of mortality models can be
chosen for the baseline mortality improvement trend. For instance, when the
focus is only on modeling the mortality for the pensioner age range, a valid
option would be to use a mortality model from the CBDX family (Dowd et al.,
2020).

##### Identifiability constraints.

To ensure model identifiability, we need to impose parameter constraints in
the model specifications of Equation (3.1). Specifically, we first impose
straightforward extensions of the identifiability constraints applied in a
classic Lee-Carter model:444Since we do not include a country-specific age-
dependent intercept, we do not impose a location constraint on the
$\kappa_{t}$-parameter.

$\displaystyle\sum_{x\>\in\>\mathcal{X}}$
$\displaystyle\left(B_{x}^{(i)}\right)^{2}=1,\hskip
14.22636pt\sum_{t\>\in\>\mathcal{T}}K_{t}^{(i)}=0,\hskip 28.45274pt$
$\displaystyle\text{for}\>i\in\\{1,\ldots,m\\}$ (3.2) $\displaystyle\hskip
28.45274pt\sum_{x\>\in\>\mathcal{X}}\left(\beta_{x}^{(j,c)}\right)^{2}=1,$
$\displaystyle\text{for}\>j\in\\{1,\ldots,l\\}\>.$

However, as we show in Appendix A, these constraints are not sufficient to
make the mortality improvement model fully identifiable. To address this
issue, we rely on additional orthogonality constraints proposed by Hunt and
Blake (2020). Specifically, for $m=2$, we impose:

$\displaystyle\displaystyle\sum_{x\>\in\>\mathcal{X}}B_{x}^{(1)}B_{x}^{(2)}=0,\hskip
14.22636pt\displaystyle\sum_{t\>\in\>\mathcal{T}}K_{t}^{(1)}K_{t}^{(2)}=0.$
(3.3)

We also impose similar constraints on the parameters that model the country-
specific deviation from the common mortality improvement rates. Appendix A
outlines the implementation details of the constraints in Equations (3.2) and
(3.3) in the Newton-Raphson algorithm that is used to fit the baseline
mortality improvement model.

##### Equivalent mortality model specification.

We translate the baseline mortality improvement model into a baseline
mortality model for the logarithm of the force of mortality $\mu_{x,t}^{(c)}$.
Appendix B shows that the model specification in Equation (3.1) is equivalent
to:

$\log\mu_{x,t}^{(c)}=\log\mu_{x,t_{\min}}^{(c)}+\left(t-t_{\min}\right)A_{x}+\displaystyle\sum_{i=1}^{m}B_{x}^{(i)}L_{t}^{(i)}+\displaystyle\sum_{j=1}^{l}\beta_{x}^{(j,c)}\lambda_{t}^{(j,c)},$
(3.4)

with $L_{t_{\min}}^{(i)}=0$ for all $i=1,...,m$,
$\lambda_{t_{\min}}^{(j,c)}=0$ for all $j=1,...,l$, and

$L_{t}^{(i)}=\displaystyle\sum_{\tau=t_{\min}+1}^{t}K_{\tau}^{(i)},\hskip
14.22636pt\lambda_{t}^{(j,c)}=\displaystyle\sum_{\tau=t_{\min}+1}^{t}\kappa_{\tau}^{(j,c)},$
(3.5)

for $t>t_{\min}$. In this expression, $\mu_{x,t_{\min}}^{(c)}$ is the force of
mortality at age $x$ in the initial year of the calibration period, typically
estimated by the crude central death rate
$\hat{m}_{x,t_{\min}}^{(c)}=d_{x,t_{\min}}^{(c)}/E_{x,t_{\min}}^{(c)}$.
Appendix B discusses the corresponding identifiability constraints of the
model specified for the (log-transformed) force of mortality in Equation
(3.4).

##### Distributional assumptions.

Hunt and Villegas (2021) present two methods for calibrating a mortality
improvement model under a Poisson assumption for the number of deaths. As
suggested by Hunt and Villegas (2021), we follow the indirect estimation
procedure for stability and reliability reasons. This indirect approach relies
on the model specification in Equation (3.4). We assume:

$\displaystyle
D_{x,t}^{(c)}\sim\text{Pois}\left(E_{x,t}^{(c)}\cdot\mu_{x,t_{\min}}^{(c)}\cdot\exp\left(\left(t-t_{\min}\right)A_{x}+\displaystyle\sum_{i=1}^{m}B_{x}^{(i)}L_{t}^{(i)}+\displaystyle\sum_{j=1}^{l}\beta_{x}^{(j,c)}\lambda_{t}^{(j,c)}\right)\right).$
(3.6)

We estimate the involved parameters using maximum likelihood (ML) estimation.
In Section 4.1, we further outline the calibration strategy using the
distributional assumption in Equation (3.6).

###  3.2 The regime-switching model with age-specific effect

##### Model specification.

To construct the baseline model’s residuals, we replace the mortality
improvement rates in Equation (2.4) with their empirical counterpart and
subtract the fitted multi-population model, as specified in Equation (3.1):

$\displaystyle
z_{x,t}^{(c)}:=\log\hat{m}_{x,t}^{(c)}-\log\hat{m}_{x,t-1}^{(c)}-\left(\hat{A}_{x}+\displaystyle\sum_{i=1}^{m}\hat{B}_{x}^{(i)}\hat{K}_{t}^{(i)}+\displaystyle\sum_{j=1}^{l}\hat{\beta}_{x}^{(j,c)}\hat{\kappa}_{t}^{(j,c)}\right),$

with $t\in\mathcal{T}\backslash\left\\{t_{\min}\right\\}$, $x\in\mathcal{X}$.
We assume that the $z_{x,t}^{(c)}$ are realizations of a random variable
$Z_{x,t}^{(c)}$, defined as:

$\displaystyle
Z_{x,t}^{(c)}:=\mathfrak{B}_{x}^{(c)}Y_{t}^{(c)}+\epsilon_{x,t}^{(c)},$ (3.7)

where the value of $Y_{t}^{(c)}$ depends on the state of an underlying Markov
chain ${\rho_{t}^{(c)}}$. At each time point $t$, $\rho_{t}^{(c)}$ can take
one of two possible values, LVS (low volatility state) or HVS (high volatility
state), which determines the value of $Y_{t}^{(c)}$. In the LVS, $Y_{t}^{(c)}$
is equal to zero, such that the mortality dynamics are described by the
baseline mortality improvement model. However, in the HVS, $Y_{t}^{(c)}$
follows a normal distribution with mean $\mu_{H}$ and variance
$\sigma_{H}^{2}$, allowing for normally distributed mortality shocks that
affect the mortality improvement trend at the country-specific level. We
define:

$\displaystyle Y_{t}^{(c)}\sim\begin{cases}0&\hskip
28.45274pt\text{if}\>\rho_{t}^{(c)}=\text{LVS}\\\
\mathcal{N}(\mu_{H},\sigma_{H}^{2})&\hskip
28.45274pt\text{if}\>\rho_{t}^{(c)}=\text{HVS}.\end{cases}$ (3.8)

Empirical evidence suggests that the impact of mortality shocks on mortality
(improvement) rates is not uniform across all age groups, as we witnessed
e.g., throughout the COVID-19 pandemic. Therefore, we model the age-specific
impact of the mortality shocks with a separate age effect
$\mathfrak{B}_{x}^{(c)}$. The country-specific residuals
$\epsilon_{x,t}^{(c)}$ represent the age- and time-specific variations in the
mortality improvement rates that are not captured by the baseline and neither
by the regime-switching model. We assume that these residuals
$\epsilon_{x,t}^{(c)}$ are independent and normally distributed with mean $0$
and standard deviation $\sigma_{e}(x,t)$.555We allow $\sigma_{e}(x,t)$ to be
dependent on age $x$ and time $t$. This assumption may be relaxed to a
constant standard deviation. In addition, we assume independence between
$Y_{t}^{(c)}$ and $\epsilon_{x,t}^{(c)}$ for all $x\in\mathcal{X}$ and
$t\in\mathcal{T}$.

The distribution of $Z_{x,t}^{(c)}$ is regime-switching and alternates as
follows:

$\displaystyle
Z_{x,t}^{(c)}\sim\begin{cases}\mathcal{N}\left(0,\sigma_{e}^{2}(x,t)\right)&\hskip
28.45274pt\text{if}\>\rho_{t}^{(c)}=\text{LVS}\\\
\mathcal{N}\left(\mathfrak{B}_{x}^{(c)}\mu_{H},\left(\mathfrak{B}_{x}^{(c)}\right)^{2}\sigma_{H}^{2}+\sigma^{2}_{e}(x,t)\right)&\hskip
28.45274pt\text{if}\>\rho_{t}^{(c)}=\text{HVS}.\end{cases}$ (3.9)

##### Identifiability constraints.

To make the model specifications of the regime-switching model identifiable,
we additionally require that:666Multiplying $\mathfrak{B}_{x}^{(c)}$ by a
constant factor $C$ and dividing both the mean and standard deviation of
$Y_{t}^{(c)}$ by the same constant $C$ yields an equivalent model
specification as described in Equation (3.7). To ensure model identifiability,
we therefore use the same type of constraint as the one imposed on the age
effects in Equation (3.2).

$\displaystyle\sum_{x\in\mathcal{X}}\left(\mathfrak{B}_{x}^{(c)}\right)^{2}=1.$

##  4 Calibration strategy

###  4.1 Calibrating the baseline mortality improvement model

In Section 4.1.1, we first calibrate the baseline mortality improvement model
on the complete calibration period $\mathcal{T}$. Then, Section 4.1.2 performs
an outlier detection strategy on the calibrated common period effects. This
allows us to identify the years in which mortality shocks took place. We
remove the observations corresponding to these years and recalibrate the
baseline mortality model. As such, we obtain the baseline improvement model
that captures the overall country-specific mortality decline.

####  4.1.1 Calibrating the baseline mortality model on the complete
calibration period

We calibrate the model in Equation (3.6) subject to the identifiability
constraints in Equations (B.3) and (B.4) through a two-step procedure. The
first step in this procedure estimates the parameters in the common trend,
while the second step puts focus on the parameters in the country-specific
deviation from that trend.

  1. 1.

In the first step, we aggregate the number of deaths and exposures by age
$x\in\mathcal{X}$ and year $t\in\mathcal{T}$ across the set $\mathcal{C}$ of
countries that constitute the common, multi-population trend. We denote the
aggregated age and period-specific deaths and exposures as:

$\displaystyle
D_{x,t}^{(A)}=\displaystyle\sum_{c\>\in\>\mathcal{C}}D_{x,t}^{(c)},\hskip
14.22636ptE_{x,t}^{(A)}=\displaystyle\sum_{c\>\in\>\mathcal{C}}E_{x,t}^{(c)}.$

We assume that the aggregated deaths random variable $D_{x,t}^{(A)}$ follows a
Poisson distribution with mean $E_{x,t}^{(A)}\cdot\exp(\eta_{x,t}^{(A)})$,
with

$\displaystyle\eta_{x,t}^{(A)}:=\log\mu_{x,t_{\min}}^{(A)}+\left(t-t_{\min}\right)A_{x}+\displaystyle\sum_{i=1}^{m}B_{x}^{(i)}L_{t}^{(i)},$

where $\mu_{x,t_{\min}}^{(A)}$ is replaced by its empirical counterpart
$\hat{m}_{x,t_{\min}}^{(A)}:=d_{x,t_{\min}}^{(A)}/E_{x,t_{\min}}^{(A)}$, i.e.
the observed aggregated common death rate at age $x$ in the initial year of
the calibration period $t_{\min}$. We estimate the common trend parameters by
maximizing the following Poisson log-likelihood:

$\displaystyle\displaystyle\max_{A_{x},\\{B_{x}^{(i)}\\}_{i},\\{L_{t}^{(i)}\\}_{i}}\>\displaystyle\sum_{t\in\mathcal{T}}\sum_{x\in\mathcal{X}}\left(d_{x,t}^{(A)}\eta_{x,t}^{(A)}-E_{x,t}^{(A)}\cdot\exp\left(\eta_{x,t}^{(A)}\right)\right),$

where we omit the terms in the log-likelihood that do not depend on the
parameters of interest. We maximize this log-likelihood using the Newton-
Raphson algorithm subject to the constraints on the common age and period
effects, as given in Equations (B.3) and (B.4).

  2. 2.

In the second step, we assume that the number of deaths random variable for
the country of interest $c$, i.e. $D_{x,t}^{(c)}$, follows a Poisson
distribution with mean $E_{x,t}^{(c)}\cdot\exp(\eta_{x,t}^{(c)})$. Since we
already calibrated the parameters in the common part of the mortality model,
we write:

$\displaystyle\eta_{x,t}^{(c)}=\hat{\eta}_{x,t}^{(A)}+\eta_{x,t}^{(\text{dev})},$
(4.1)

with

$\displaystyle\eta_{x,t}^{(\text{dev})}:=\log\left(\frac{\mu_{x,t_{\min}}^{(c)}}{\mu_{x,t_{\min}}^{(A)}}\right)+\displaystyle\sum_{j=1}^{l}\beta_{x}^{(j,c)}\lambda_{t}^{(j,c)},$
(4.2)

where $\mu_{x,t_{\min}}^{(c)}$ and $\mu_{x,t_{\min}}^{(A)}$ are replaced by
their empirical counterparts. We then estimate the country-specific parameters
by maximizing the following Poisson log-likelihood:

$\displaystyle\displaystyle\max_{\\{\beta_{x}^{(j,c)}\\}_{j},\\{\lambda_{t}^{(j,c)}\\}_{j}}\>\displaystyle\sum_{t\in\mathcal{T}}\sum_{x\in\mathcal{X}}\left(d_{x,t}^{(c)}\eta_{x,t}^{(c)}-E_{x,t}^{(c)}\cdot\exp\left(\eta_{x,t}^{(c)}\right)\right),$

subject to the constraints on the country-specific age and period effects in
Equations (B.3) and (B.4).

####  4.1.2 Towards an outlier-free baseline mortality improvement model

The quality of the fit of the baseline mortality improvement model may be
influenced significantly by the presence of mortality shocks in the mortality
data. Hereto, we investigate the presence of outliers in the calibrated common
period effects $\hat{L}_{t}^{(i)}$, where $i\in\\{1,\ldots,m\\}$. Our focus on
outlier detection in the common trend parameters only is motivated in the case
study in Section 6.

Consider a calibrated period effect $\hat{L}_{t}^{(i)}$ and denote it as
$L_{t}^{(i)}$ for notational convenience. We decompose the time series
$L_{t}^{(i)}$ as (Cleveland and Tiao, 1976):

$L_{t}^{(i)}=T_{t}^{(i)}+R_{t}^{(i)},$

for $i\in\\{1,2,...,m\\}$, where $T_{t}^{(i)}$ is a trend component and
$R_{t}^{(i)}$ a so-called remainder component.777Cleveland and Tiao (1976)
also consider a seasonal component. Since we work with annual data, we ignore
the seasonal component. After estimating the trend component, the remainder
component equals:

$\hat{R}_{t}^{(i)}=L_{t}^{(i)}-\hat{T}^{(i)}_{t}.$

Using an $m$-dimensional outlier detection technique, we determine the
outliers in the set of remainder components:

$\mathcal{R}=\left\\{\hat{\boldsymbol{R}}_{t}:=\left(\hat{R}_{t}^{(1)},\hat{R}_{t}^{(2)},\ldots,\hat{R}_{t}^{(m)}\right)\mid
t\in\mathcal{T}\right\\}.$

Appendix C explains our proposed strategy to detect the outlying years in
$\mathcal{R}$ using the Mahalanobis distance. We denote the set of years that
are identified as outlying with $\mathcal{T}^{o}$.

Next, we recalibrate the baseline mortality model in Equation (3.6) on the
outlier-free calibration period
$\widetilde{\mathcal{T}}=\mathcal{T}\backslash\left\\{\mathcal{T}^{o}\right\\}$.
As a result, the recalibrated period effects $\hat{L}_{t}^{(i)}$ will be
missing for the years detected as outlying. We account for this using missing
value imputation.888We explain the details of the missing value imputation
technique in Section 6.2.2 of the case study. We denote the final recalibrated
parameters in our baseline mortality model as $\hat{A}_{x}$,
$\hat{B}_{x}^{(1)}$, …, $\hat{B}_{x}^{(m)}$, $\hat{L}_{t}^{(1)}$, …,
$\hat{L}_{t}^{(m)}$, $\hat{\beta}_{x}^{(1,c)}$, …, $\hat{\beta}_{x}^{(l,c)}$,
$\hat{\lambda}_{t}^{(1,c)}$, …, $\hat{\lambda}_{t}^{(l,c)}$.

After the recalibration of the age and period effects in the baseline
mortality model, we recover the period effects of the original mortality
improvement model via back-transformation starting from Equation (3.5):

$\displaystyle\hat{K}_{t}^{(i)}=\hat{L}_{t}^{(i)}-\hat{L}_{t-1}^{(i)},\hskip
14.22636pt\hat{\kappa}_{t}^{(j,c)}=\hat{\lambda}_{t}^{(j,c)}-\lambda_{t-1}^{(j,c)},$
(4.3)

for $t\in\mathcal{T}$, $i\in\\{1,..,m\\}$ and $j\in\\{1,...,l\\}$. The age
effects in the original model specifications remain unchanged. Note that, by
construction, the identifiability constraints in Equations (B.3) and (B.4) are
fulfilled. In this way, a baseline mortality improvement model that captures
the overall country-specific mortality decline results. We denote the final
calibrated baseline improvement model as:

$\displaystyle\log\hat{\mu}_{x,t}^{(c)}-\log\hat{\mu}_{x,t-1}^{(c)}=\hat{A}_{x}+\displaystyle\sum_{i=1}^{m}\hat{B}_{x}^{(i)}\hat{K}_{t}^{(i)}+\displaystyle\sum_{j=1}^{l}\hat{\beta}_{x}^{(j,c)}\hat{\kappa}_{t}^{(j,c)}.$
(4.4)

The trend line in Figure 1 visualizes the decrement in mortality rates
attributed to the fitted baseline improvement model in Equation (4.4).

###  4.2 Calibrating the regime-switching model for age-specific mortality
shocks

We now explain the calibration strategy for the regime-switching model that
drives the mortality shocks and their age-specific effects, see Section 3.2.
In the low volatility regime, the dynamics of mortality are modeled by the
baseline mortality improvement model, which was calibrated in Section 4.1. In
the high volatility regime, the mortality trend experiences disruptions
through normally distributed and age-specific mortality shocks.

Starting from the baseline model’s residuals $Z_{x,t}$, as defined in Equation
(3.7), we switch to vector notation:

$\displaystyle\boldsymbol{Z}_{t}^{(c)}=\boldsymbol{\mathfrak{B}}^{(c)}Y_{t}^{(c)}+\boldsymbol{E}_{t}^{(c)},$

where $\boldsymbol{Z}_{t}^{(c)}:=(Z_{x,t}^{(c)})_{x\in\mathcal{X}}$. We denote
$\boldsymbol{\mathfrak{B}}^{(c)}:=(\mathfrak{B}_{x}^{(c)})_{x\in\mathcal{X}}$
for the vector of the mortality shocks’ age-specific effects and
$\boldsymbol{E}_{t}^{(c)}=(\epsilon_{x,t}^{(c)})_{x\in\mathcal{X}}$ for the
vector of error terms. We obtain:999$|\mathcal{X}|$ denotes the number of ages
in the age range $\mathcal{X}$.

$\displaystyle\boldsymbol{Z}_{t}^{(c)}$
$\displaystyle\sim\begin{cases}\mathcal{N}_{|\mathcal{X}|}\left(\boldsymbol{0},\sigma_{e}^{2}(x,t)\boldsymbol{I}_{|\mathcal{X}|}\right)\hskip
28.45274pt&\text{if}\>\rho_{t}^{(c)}=\text{LVS}\\\
\mathcal{N}_{|\mathcal{X}|}\left(\boldsymbol{\mathfrak{B}}^{(c)}\mu_{H},\boldsymbol{\mathfrak{B}}^{(c)}\left(\boldsymbol{\mathfrak{B}}^{(c)}\right)^{T}\sigma_{H}^{2}+\sigma_{e}^{2}(x,t)\boldsymbol{I}_{|\mathcal{X}|}\right)\hskip
28.45274pt&\text{if}\>\rho_{t}^{(c)}=\text{HVS},\end{cases}$ (4.5)

with $\boldsymbol{I}_{|\mathcal{X}|}$ the identity matrix of size
$|\mathcal{X}|$ and $\mathcal{N}_{|\mathcal{X}|}(\cdot,\cdot)$ the
$|\mathcal{X}|$-dimensional Gaussian distribution. We denote the time-
invariant probabilities to transition between the states of the Markov chain
$(\rho_{t}^{(c)})_{t\in\mathcal{T}}$ as:

$\displaystyle
p_{ij}=\mathbb{P}\left(\rho_{t}^{(c)}=j\mid\rho_{t-1}^{(c)}=i\right),\hskip
14.22636pti,j\in\\{1,2\\}\>\>\text{and}\>\>t\in\mathcal{T}\backslash\\{t_{\min}\\},$
(4.6)

where $p_{11}+p_{12}=1$ and $p_{21}+p_{22}=1$.101010For notational
convenience, state ‘1’ refers to the LVS and state ‘2’ refers to the HVS in
the sequel of this paper. If a mortality shock occurs in a given year $t$, we
typically observe a positive outlier in the mortality improvement rates in
year $t$, i.e. for $\log\hat{\mu}_{x,t}^{(c)}-\log\hat{\mu}_{x,t-1}^{(c)}$,
and a negative outlier in year $t+1$, i.e. for
$\log\hat{\mu}_{x,t+1}^{(c)}-\log\hat{\mu}_{x,t}^{(c)}$. Hence, at least two
consecutive observations of the mortality improvement rates are impacted by
the shock. To accommodate this observation, we introduce a memory state into
the Markov chain, as detailed in Appendix E.1. As such, we guarantee that the
Markov chain remains in the HVS for a minimum duration of two years once it
transitions to that state. In this way, we better mimic situations that occur
in real-world settings.

The parameter vector in the regime-switching model is $|\mathcal{X}|+5$
dimensional and will be denoted as:111111By default we assume a constant
standard deviation $\sigma_{e}$ for the residuals $\epsilon_{x,t}$ in the
mortality improvement model with shock regime. However, in Section 6, we relax
this assumption and tailor it to the data at hand.

$\displaystyle\Theta=\left\\{\mu_{H},\sigma_{H},\sigma_{e},p_{12},p_{21},\boldsymbol{\mathfrak{B}}^{(c)}\right\\}\in\mathbb{R}^{|\mathcal{X}|+5}.$

##### Constructing the log-likelihood.

We obtain the observed sample $\\{\boldsymbol{z}_{t}\\}_{t\in\mathcal{T}}$ as
the residuals from the calibrated baseline mortality improvement model in
Equation (4.4). The log-likelihood of the sample equals:

$\displaystyle l(\Theta)$ $\displaystyle=\log
f\left(\boldsymbol{z}_{t_{\min}}^{(c)},\boldsymbol{z}_{t_{\min+1}}^{(c)},...,\boldsymbol{z}_{t_{\max}}^{(c)}\mid\Theta\right)$
(4.7) $\displaystyle=\sum_{t\in\mathcal{T}}\log
f\left(\boldsymbol{z}_{t}^{(c)}\mid\boldsymbol{z}_{t-1}^{(c)},...,\boldsymbol{z}_{t_{\min}}^{(c)},\Theta\right).$

Appendix E demonstrates the calculation of the involved conditional
probabilities in a recursive way. Additionally, we explain our strategy to
determine starting values for the recursion.

##### Optimizing the log-likelihood.

Using the starting values and log-likelihood derived in Appendix E, we obtain
the optimal set of parameters as:

$\hat{\Theta}=\operatorname*{arg\,max}_{{\Theta}}\>l(\Theta).$

##  5 Mortality projection strategy

Next, we proceed with the mortality projection strategy. We first establish
the time series models for the calibrated common period effects
($\hat{K}_{t}^{(1)}$, …, $\hat{K}_{t}^{(m)}$) and the country-specific period
effects ($\hat{\kappa}_{t}^{(1,c)}$, …, $\hat{\kappa}_{t}^{(l,c)}$) for
$t\in\mathcal{T}$. We use a multivariate time series model for the joint
projection of both the common and country-specific period effects over time.
Then, we generate the states that the Markov chain $\rho_{t}^{(c)}$ will
occupy in the future as well as the severity $Y_{t}^{(c)}$ of future mortality
shocks in the high volatility periods of this Markov chain. We combine the
projected time series and mortality shocks with age-specific effects (see
Equations (3.1) and (3.7)) and obtain as such the projections for the future
mortality rates.

###  5.1 Specifying the time series for the calibrated period effects

We represent the vector of period effects as:

$\displaystyle\boldsymbol{K}_{t}=\begin{pmatrix}K_{t}^{(1)}\\\ \vdots\\\
K_{t}^{(m)}\end{pmatrix}\in\mathbb{R}^{m\times 1},\hskip
14.22636pt\boldsymbol{\kappa}_{t}^{(c)}=\begin{pmatrix}\kappa_{t}^{(1,c)}\\\
\vdots\\\ \kappa_{t}^{(l,c)}\end{pmatrix}\in\mathbb{R}^{l\times 1},$

and are interested in a multivariate time series model for these combined
period effects. We denote:

$\displaystyle\boldsymbol{\mathcal{K}}_{t}=\begin{pmatrix}\boldsymbol{K}_{t}\\\
\boldsymbol{\kappa}_{t}^{(c)}\end{pmatrix}\in\mathbb{R}^{(m+l)\times 1}.$

Given our focus on a mortality improvement model, as motivated in Section 1,
we impose a VARMA($p$,$q$) model on $\boldsymbol{\mathcal{K}}_{t}$:121212Due
to our focus on mortality improvement rates, the period effects have already
been implicitly differenced, as demonstrated in Equation (4.3).

$\displaystyle\boldsymbol{\mathcal{K}}_{t}$
$\displaystyle=\boldsymbol{c}+\sum_{r=1}^{p}\boldsymbol{\Phi}_{r}\boldsymbol{\mathcal{K}}_{t-r}+\sum_{s=1}^{q}\boldsymbol{\Psi}_{s}\boldsymbol{\delta}_{t-s}+\boldsymbol{\delta}_{t},$
(5.1) $\displaystyle\boldsymbol{\delta}_{t}$
$\displaystyle\sim\mathcal{N}_{m+l}\left(\boldsymbol{0},\boldsymbol{\Sigma}_{\delta}\right),$

for any $t\in\mathcal{T}$ and where $\boldsymbol{\Phi}_{r}$ and
$\boldsymbol{\Psi}_{s}$ for all $r,s$ are $(m+l)$-dimensional matrices.
Additionally, $\boldsymbol{\delta}_{t}$ denotes the joint $(m+l)$-dimensional
vector of error terms at time $t$. We estimate the parameters in Equation
(5.1) using maximum likelihood estimation. The case study in Section 6 adopts
a weighting scheme in the maximum likelihood estimation of the VARMA model. As
such, the modeler can put more weight on recent observations to better reflect
the recent trends in the mortality decline.

###  5.2 Mortality rate projections

We are now ready to generate a single trajectory for the mortality rates,
$q_{x,t}^{(c)}$, for each age $x$ over the projection period
$t\in\mathcal{T}^{\text{{\tiny pred}}}:=\\{t_{\max}+1,\ldots,T\\}$. We repeat
this procedure $10\ 000$ times to assess the uncertainty in the projections.
We denote each trajectory as $\hat{q}_{x,t,\iota}^{(c)}$, with
$x\in\mathcal{X}$, $t\in\mathcal{T}^{\text{{\tiny pred}}}$ and $\iota$ the
number of the generated path. Hereto, we first generate trajectories for both
the estimated time series models from Section 5.1 as well as the Markov chain
regime-switching model with estimated transition probabilities $\hat{p}_{ij}$
from Section 4.2. We proceed as follows:

  1. 1.

By generating new error terms $\boldsymbol{\delta}_{t,\iota}$ for
$t\in\mathcal{T}^{\text{{\tiny pred}}}$ from the fitted multivariate normal
distribution with mean $\boldsymbol{0}$ and covariance matrix
$\hat{\boldsymbol{\Sigma}}_{\delta}$, we obtain a trajectory for the future
$\hat{\boldsymbol{\mathcal{K}}}_{t,\iota}$ using the calibrated time series
model from Equation (5.1). This leads to generated trajectories for the common
period effects $\hat{K}_{t,\iota}^{(i)}$ for $i=1,..,,m$ and the country-
specific period effects $\hat{\kappa}_{t,\iota}^{(j,c)}$ for $j=1,...,l$.

  2. 2.

To generate a trajectory for the regime-switching model over the projection
period $\mathcal{T}^{\text{{\tiny pred}}}$, we first generate a trajectory for
the underlying Markov chain $\rho_{t}^{(c)}\in\\{1,2\\}$ for
$t\in\mathcal{T}^{\text{{\tiny pred}}}$. We initialize the trajectory of the
Markov chain at time $t_{\max}$ with the most probable state as derived from
the calibrated regime-switching model.131313This is the state $r$ with the
largest estimated probability
$\mathbb{P}(\rho_{t_{\max}}^{(c)}=r\mid\boldsymbol{z}_{t_{\max}}^{(c)},\ldots,\boldsymbol{z}_{t_{\min}}^{(c)},\hat{\boldsymbol{\Theta}})$.
Based on this initial state and the calibrated transition probabilities, we
can then generate the state the Markov chain will occupy at any future time
point $t\in\mathcal{T}^{\text{{\tiny pred}}}$. This enables us to generate a
trajectory for $\hat{Y}_{t,\iota}^{(c)}$ using the fitted normal distribution
in the high volatility regime with mean $0$ and variance
$\hat{\sigma}_{H}^{2}$, see Equation (3.8). Furthermore, we ensure that the
generated mortality shocks lead to sudden increases in the mortality rates,
possibly spanning multiple years, followed by a return to the overall baseline
mortality decline, see Appendix E.2 for more details.

We generate 10 000 trajectories for the force of mortality of a person aged
$x\in\mathcal{X}$ over projection period $t\in\mathcal{T}^{\text{{\tiny
pred}}}$ in a recursive way:

$\displaystyle\log\hat{\mu}_{x,t,\iota}^{(c)}=\log\hat{\mu}_{x,t-1,\iota}^{(c)}+\displaystyle\sum_{i=1}^{m}\hat{B}_{x}^{(i)}\hat{K}_{t,\iota}^{(i)}+\displaystyle\sum_{j=1}^{l}\hat{\beta}_{x}^{(j,c)}\hat{\kappa}_{t,\iota}^{(j,c)}+\hat{\mathfrak{B}}_{x}^{(c)}\hat{Y}_{t,\iota}^{(c)},$
(5.2)

where $\iota\in\\{1,2,\dots,10\ 000\\}$, the generated paths. We start the
recursion with the observed crude central death rates at time $t_{\max}$, i.e.
$\hat{\mu}_{x,t_{\max},\iota}^{(c)}=d_{x,t_{\max}}^{(c)}/E_{x,t_{\max}}^{(c)}$
for all $\iota$. The $\iota$-th generated mortality rate at age $x$ and time
$t$ then equals:

$\displaystyle\hat{q}_{x,t,\iota}^{(c)}$
$\displaystyle=1-\exp\left(-\hat{\mu}_{x,t,\iota}^{(c)}\right).$

##  6 Case study

###  6.1 Multi-population mortality data set

We use a multi-population mortality data set composed of a pre-selected set of
countries. Hereto we rely on the set of countries selected by the Institute of
Actuaries in Belgium (Antonio et al., 2020) and the Dutch Actuarial
Association (Royal Dutch Actuarial Association, 2020, 2022): Austria (AT),
Belgium (BE), Denmark (DK), Germany (DE), Finland (FI), France (FR), Ireland
(IE), Iceland (IS), Luxembourg (LU), the Netherlands (NL), Norway (NO), Sweden
(SE), Switzerland (CH), and the United Kingdom (UK).141414Ireland is excluded
from our case study due to the unavailability of weekly mortality information,
as elaborated later on in this section.

Whereas Antonio et al. (2020) let the calibration period of their proposed Li-
Lee model start in 1988 and Royal Dutch Actuarial Association (2020, 2022) in
1970, this study puts focus on a longer calibration period that starts in the
year 1850. This enables access to more historical mortality shocks. However,
it is not possible to retrieve mortality statistics for all thirteen
considered European countries from the year 1850 onwards. To resolve this
issue we work with a time-varying composition of the multi-population data set
and sequentially add countries when more information becomes available. Figure
4 colors the thirteen considered European countries and indicates the year
from which each country is included in the multi-population mortality trend.
We denote the set of countries that compose this mortality trend at year $t$
as $\mathcal{C}_{t}$. As an example, $\mathcal{C}_{1900}$ consists of nine
countries: AT, BE, DK, FI, FR, IS, NL, NO and SE.

Figure 4: Visualization of the thirteen European countries that are part of
the common mortality improvement trend and indicate the year when mortality
information becomes available..

As such, our resulting calibration period is
$\mathcal{T}=\\{1850,1851,...,2021\\}$ and we focus on the age range
$\mathcal{X}=\\{20,21,...,85\\}$. The exclusion of the oldest ages is a result
of the limited exposure available for these ages and is motivated by Zhou et
al. (2013) and Liu and Li (2015). Furthermore, we exclude the youngest ages
due to the high variability and limited number of deaths in this age range
which can potentially result in unstable outcomes, as discussed in Zhou et al.
(2013) and Hainaut (2012). In this case study, our focus is on modeling the
Dutch male mortality rates.151515We analyze male mortality data due to the
historical prevalence of more male victims in world wars, making the male
mortality data set particularly relevant for this paper. Of course, a similar
analysis can be performed on female or unisex data.

We retrieve the annual deaths, $d_{x,t}^{(c)}$, and exposures,
$E_{x,t}^{(c)}$, at individual ages $x\in\mathcal{X}$ for each country $c$
from the Human Mortality Database (HMD).161616The Human Mortality Database
(HMD) can be accessed at https://www.mortality.org/. We further complement the
mortality data with annual mortality statistics from Eurostat when
possible.171717The Eurostat mortality data-base can be consulted on
https://ec.europa.eu/eurostat. Following the approach from Robben et al.
(2022), we complement the data sets, when necessary, with deaths and exposures
up to the year 2021 using weekly mortality information from Eurostat and the
Short Term Mortality Fluctuation (STMF) data series on the HMD.

###  6.2 Calibrating the baseline mortality improvement model

Using the strategy outlined in Section 4.1, we calibrate the baseline
mortality improvement model over the complete calibration period
$\mathcal{T}$. In this case study, we specify a two-factor age/period
structure for both the common and the country-specific deviation trend, i.e.
we use $m=l=2$ in the baseline mortality improvement model in Equation
(3.1):181818We omit the superscript $c$ in the country-specific parameters for
notational purposes.

$\displaystyle\log\mu_{x,t}^{({\scriptscriptstyle\text{NL}})}-\log\mu_{x,t-1}^{({\scriptscriptstyle\text{NL}})}$
$\displaystyle=A_{x}+B_{x}^{(1)}K_{t}^{(1)}+B_{x}^{(2)}K_{t}^{(2)}+\beta_{x}^{(1)}\kappa_{t}^{(1)}+B_{x}^{(2)}\kappa_{t}^{(2)},$
(6.1)

and in the equivalent baseline mortality model in Equation (3.4):

$\displaystyle\log\mu_{x,t}^{({\scriptscriptstyle\text{NL}})}$
$\displaystyle=\log\mu_{x,t_{\min}}^{({\scriptscriptstyle\text{NL}})}+\left(t-t_{\min}\right)A_{x}+B_{x}^{(1)}L_{t}^{(1)}+B_{x}^{(2)}L_{t}^{(2)}+\beta_{x}^{(1)}\lambda_{t}^{(1)}+\beta_{x}^{(2)}\lambda_{t}^{(2)},$
(6.2)

where $L_{t}^{(1)}$, $L_{t}^{(2)}$, $\lambda_{t}^{(1)}$ and
$\lambda_{t}^{(2)}$ are defined as in Equation (3.5).

####  6.2.1 Detecting outliers

Figure 5 displays the calibrated parameters in the common part of the baseline
mortality model. We observe clear outliers in the calibrated common period
effects $\hat{L}_{t}^{(1)}$ and $\hat{L}_{t}^{(2)}$, representing mortality
shocks.191919For the remaining of this paper, we will use the terms ‘outliers’
and ‘mortality shocks’ interchangeably. We observe the most pronounced
outliers during the two world wars.

Figure 5: The calibrated age (top) and period (bottom) effects in the common
trend of the baseline mortality model in Equation (3.4). Calibration period
$1850$-$2021$, age range 20-85, male data.

Following Section 4.1.2, we first determine the trend component in each of the
calibrated common period effects using robust cubic smoothing splines.202020We
use the function robustSmoothPline from the aroma.light package in R, which is
a robust version of the commonly used smooth.spline function from the stats
package. Specifically, for $i=1,2$ we aim to find a function $f_{i}$ that
minimizes the following objective function:

$\displaystyle\displaystyle\hat{f}_{i}^{\lambda}=\operatorname*{arg\,min}_{f_{i}}\displaystyle\sum_{t=1850}^{2021}\left(L_{t}^{(i)}-f_{i}(t)\right)^{2}+\lambda\displaystyle\int\left(f_{i}^{\prime\prime}(t)\right)^{2}dt.$
(6.3)

We approximate the function $f_{i}(\cdot)$ with a cubic smoothing spline with
the years in the calibration period $\mathcal{T}$ as knots. To ensure the
smoothness of the function $f_{i}$, we employ a penalty function based on the
second-order derivative of $f_{i}$. The smoothing parameter $\lambda$ balances
the goodness of the spline fit and the degree of smoothness of the cubic
spline. We select $\lambda$ via the Generalized Cross Validation (GCV)
criterion (Golub et al., 1979). Figure 6 shows the resulting smoothing
splines, as obtained with the optimally selected smoothing parameter
$\hat{\lambda}_{\text{opt}}^{(i)}$.212121We a priori remove the years of the
major mortality shocks (see Figure 2) to attain better robustness of our
spline fits.

Figure 6: The robust smoothing spline fitted to the calibrated common period
effect $\hat{L}_{t}^{(1)}$ (left) and $\hat{L}_{t}^{(2)}$ (right) where we
select the optimal smoothing parameter using the GCV criterion. Calibration
period $1850$-$2021$, age range 20-85, male data.

Next, for each $i=1,2$, we define the remainder component of the calibrated
period effect $\hat{L}_{t}^{(i)}$ after removing the trend captured by the
smoothing spline:

$R_{t}^{(i)}=\hat{L}_{t}^{(i)}-\hat{f}_{i}^{\hat{\lambda}_{\text{opt}}^{(i)}}.$

We consider the bivariate time series of these remainder components:

$\mathcal{R}=\left\\{\boldsymbol{R}_{t}:=\left(R_{t}^{(1)},R_{t}^{(2)}\right)\mid
t=1850,...,2021\right\\}.$

On these remainder components, we apply the robust bivariate outlier detection
procedure explained in Appendix C. Since the mortality rates before the year
1970 are more volatile compared to the mortality rates after 1970, we perform
separate outlier detection on the set of remainder components before and after
1970. Figure 7 shows the results. The set of detected outlying years equals:

$\mathcal{T}^{o}=\left\\{1854{\text{-}}1856,\>1859,\>1866,\>1870{\text{-}}1871,\>1914{\text{-}}1920,\>1929,\>1940{\text{-}}1945,1991,\>2020{\text{-}}2021\right\\}.$

The detected outliers demonstrate a strong correspondence with the historical
mortality shocks observed in Europe, as shown in Figures 1 and 2.

Figure 7: Outlier detection using Mahalanobis distance on the remainder
components before and after 1970. Calibration period $1850$-$2021$, age range
20-85, male data.

####  6.2.2 Constructing the baseline mortality improvement model

We now recalibrate the baseline mortality model from Equation (3.6) on the
outlier-free calibration period
$\widetilde{\mathcal{T}}:=\mathcal{T}\backslash\left\\{\mathcal{T}^{o}\right\\}$.
This ensures that the calibrated baseline mortality model is not influenced by
the mortality shocks that have been observed and detected in Section 6.2.1.
However, this procedure results in missing observations in the calibrated
period effects $\hat{L}_{t}^{(j)}$ due to the exclusion of certain years from
the calibration period. To address this issue, we intervene in the Newton-
Raphson (NR) algorithm, used to calibrate the baseline mortality model as
outlined in Section 4.1.1. More specifically, at the end of each iteration in
the NR algorithm, we perform missing value imputation for the missing years in
the period effects using an exponential weighted moving average
technique.222222We achieve this by using the na_ma function of the ImputeTS
package (Moritz and Bartz-Beielstein, 2017). This results in calibrated period
effects that are defined over the complete calibration period $\mathcal{T}$
and that fulfill the identifiability constraints.

Following Equation (4.3), we transform the fitted parameters
$\hat{L}_{t}^{(i)}$ for $i=1,2$ and $\lambda_{t}^{(j)}$ for $j=1,2$ to the
parameters $\hat{K}_{t}^{(i)}$ for $i=1,2$ and $\hat{\kappa}_{t}^{(j)}$ for
$j=1,2$ in the original baseline mortality improvement model outlined in
Equation (4.4). Figure 8 shows the age and period effects in the common trend
specification (top) and the specification of the Dutch deviation from the
European trend (bottom). In line with our expectations, the calibrated period
effects no longer exhibit any severe outliers.

Figure 8: The calibrated common (top) and country-specific (bottom) age and
period effects in the baseline mortality improvement model. Outlier-free
calibration period
$\widetilde{\mathcal{T}}=\mathcal{T}\backslash\left\\{\mathcal{T}^{o}\right\\}$
combined with missing value imputation, age range 20-85, male data.

We denote the calibrated age/period structure in the Dutch baseline mortality
improvement model (see Equation (6.1)) as:

$\displaystyle\hat{\eta}_{x,t}^{({\scriptscriptstyle\text{NL}})}:=\hat{A}_{x}+\hat{B}_{x}^{(1)}\hat{K}_{t}^{(1)}+\hat{B}_{x}^{(2)}\hat{K}_{t}^{(2)}+\hat{\beta}_{x}^{(1)}\hat{\kappa}_{t}^{(1)}+\hat{B}_{x}^{(2)}\hat{\kappa}_{t}^{(2)}.$

Following Section 3.2, we calculate the baseline model’s residuals as:

$\displaystyle
z_{x,t}:=\log\hat{m}_{x,t}^{({\scriptscriptstyle\text{NL}})}-\log\hat{m}_{x,t-1}^{({\scriptscriptstyle\text{NL}})}-\hat{\eta}_{x,t}^{({\scriptscriptstyle\text{NL}})},$
(6.4)

which are perturbed by periods of high volatility caused by mortality shocks.
Figure 9 shows the time series of observed $z_{x,t}$ for each age
$x\in\mathcal{X}$ and $t\in\mathcal{T}$.

Figure 9: The time series of realized residuals $z_{x,t}$ of the calibrated
baseline mortality improvement model for all ages $x\in\mathcal{X}$.
Calibration period 1851-2021, age range 20-85, male data.

###  6.3 Calibrating the regime-switching model

####  6.3.1 Regime-switching model with age-specific effect

We calibrate the regime-switching model separately on two distinct age groups:
$\mathcal{X}^{(1)}$, covering ages 20-59, and $\mathcal{X}^{(2)}$, covering
ages 60-85. This split enables us to distinguish the occurrence of different
types of mortality shocks. By doing so, we can effectively distinguish the
impact of events like the COVID-19 pandemic, which mainly affects older age
groups, from pandemics and wars which primarily affect younger age groups.
Therefore, we impose the following adjusted specifications for the baseline
model’s residuals:

$\displaystyle
z_{x,t}\sim\begin{cases}\mathfrak{B}_{x}^{(1)}Y_{t}^{(1)}+\epsilon_{x,t}^{(1)}&\text{for
}x\in\mathcal{X}^{(1)}\\\
\mathfrak{B}_{x}^{(2)}Y_{t}^{(2)}+\epsilon_{x,t}^{(2)}&\text{for
}x\in\mathcal{X}^{(2)},\end{cases}$ (6.5)

where $Y_{t}^{(1)}$ and $Y_{t}^{(2)}$ are both regime-switching and alternate
between a low and high volatility state, characterized by two independent
Markov chains $\rho_{t}^{(1)}$ and $\rho_{t}^{(2)}$, respectively. These
regime-switching models and Markov chains have their own set of parameters and
transition probabilities. The time series vector
$\boldsymbol{z}_{t}:=(z_{x,t})_{x\in\mathcal{X}^{(j)}}$ follows for each age
group $\mathcal{X}^{(j)}$ an alternating multivariate normal distribution, as
derived in Equation (4.5) of Section 4.2. In the remainder of this section, we
omit the superscript related to the age range under consideration for
notational purposes.

In accordance with the methodology outlined in Section 4.2, we adopt a
flexible specification for the standard deviation of the residuals
$\epsilon_{x,t}$. From Figure 1 we learn that the death rates in the years
after 1970 exhibit lower volatility than those observed in previous years.
Moreover, Figure 9 reveals that the standard deviation of the error term
$\epsilon_{x,t}$ decreases with age.232323We ignore the outliers in Figure 9
to assess the variability of $\epsilon_{x,t}$. To accommodate these
observations, we propose the following specification for $\sigma_{e}(x,t)$:

$\displaystyle\sigma_{e}(x,t)$
$\displaystyle=\sigma_{e_{1}}(x)\mathbbm{1}_{\left\\{t<1970\right\\}}+\sigma_{e_{2}}(x)\mathbbm{1}_{\left\\{t\geq
1970\right\\}},$ (6.6) $\displaystyle\sigma_{e_{j}}(x)$
$\displaystyle=\sigma_{e_{j}}+\text{slope}_{j}\cdot(x-x_{\min}),$

for $j\in\\{1,2\\}$. Further, $x_{\min}$ refers to the minimum age in the age
group under consideration and, e.g., $\mathbbm{1}_{\left\\{t<1970\right\\}}$
is one for $t$ smaller than the year 1970 and zero elsewhere. This introduces
extra flexibility regarding the variability of the residuals $\epsilon_{x,t}$
in the period before and after the year 1970. The parameter ‘slope’ captures
the linear decrease in the variability with age, before and after 1970. We use
the same specification for $\sigma_{e}(x,t)$ for both age groups
$\mathcal{X}^{(1)}$ and $\mathcal{X}^{(2)}$.

####  6.3.2 Calibration results for the regime-switching model

We provide all technical details about the calibration of the regime-switching
model in Appendix E.1. Table 1 and Figure 10 display the estimated parameters
in the regime-switching model, which we calibrate separately for the age
groups 20-59 and 60-85. We interpret the results for the age group 20-59. The
estimated value for $p_{1,2}$ indicates that there is a $4.7\%$ probability of
transitioning from the low to the high volatility state. The standard
deviation of the shock random variable $Y_{t}^{(1)}$ in the high volatility
regime, i.e. $\sigma_{H}=2.41$, is noticeably higher compared to that in the
low volatility regime.242424To compute the standard deviation in the low
volatility state, we use Equation (6.6). The calibrated age dependency
$\hat{\mathfrak{B}}_{x}$ of mortality shocks for ages 20-59 exhibits a
decreasing trend, resulting in less severe shocks for older individuals
compared to younger ones. This can be attributed to the mortality shocks
resulting from the two world wars where the majority of the victims during the
world wars were young men.  Parameter 20-59 60-85 $p_{1,2}$ 0.04709 0.06656
$p_{2,1}$ 0.34207 0.68966 $\sigma_{e_{1}}$ 0.12750 0.05823 $\text{slope}_{1}$
-0.00156 0.00069 $\sigma_{e_{2}}$ 0.18173 0.03944 $\text{slope}_{2}$ -0.00342
-0.00048 $\mu_{H}$ -0.01602 0.03155 $\sigma_{H}$ 2.40648 0.89818 x Table 1:
Calibrated parameters in the regime-switching model for the two age groups.
Calibration period 1850-2021, male data. Figure 10: Calibrated age-specific
effect $\mathfrak{B}_{x}$ in the regime-switching model for age groups 20-59
and 60-85. Calibration period 1850-2021, male data. The green line in Figure
11 shows the conditional probability of being in the low volatility state
(LVS, state 1) in year $t$, as derived in Equation (E.3). The figure reveals
alternating periods of high and low volatility. The blue dots in this figure
represent the age-averaged realized residuals $z_{x,t}$ in the calibrated
baseline mortality improvement model. The outlying dots correspond to
historical mortality shocks and are captured by the HVS. The results indicate,
for instance, that during the COVID-19 pandemic period 2020-2021 the Markov
chain occupies the HVS but only for the age group 60-85. This is in line with
our understanding that the pandemic primarily affects older individuals.

Figure 11: The green line represents the conditional probability of being in
the LVS (state 1) at year $t$, calculated using Equation (E.3). These
probabilities are calculated with the optimized parameter values in Table 1
and Figure 10. The blue dots represent the age-averaged realized residuals
$z_{x,t}$ versus year $t$. The left panel shows the results for the age group
20-59 and the right panel for the age group 60-85. Calibration period
1850-2021, male data.

###  6.4 Time series models and resulting mortality projections

We follow the projection strategy outlined in Section 5 to project mortality
rates with our proposed mortality improvement model with a shock regime.

####  6.4.1 Projecting common and country-specific period effects using a
multivariate time series model

In accordance with Section 5.1 and inspired by the random behavior of the
calibrated period effects illustrated in Figure 8, we propose the following
multivariate time series model for the four-dimensional vector of period
effects
$\boldsymbol{\mathcal{K}}_{t}=\left(K_{t}^{(1)},K_{t}^{(2)},\kappa_{t}^{(1)},\kappa_{t}^{(2)}\right)$
for $t\in\\{1851,1852,\ldots,2021\\}$:

$\displaystyle\boldsymbol{\mathcal{K}}_{t}=\boldsymbol{c}+\boldsymbol{W}_{t},$
(6.7)

where $\boldsymbol{c}=[c_{1},c_{2},0,0]^{T}$ are the mean parameters or
intercepts and $\boldsymbol{W}_{t}$ is a four-dimensional vector of white
noise terms with mean $\boldsymbol{0}$ and covariance matrix
$\boldsymbol{\Sigma}_{w}$.252525We use zero intercepts in the time series
models of the country-specific period effects to avoid divergence from the
common improvement trend. To estimate the time series models in Equation
(6.7), we use a weighted variant of maximum likelihood estimation, see
Appendix D for more details. This weighting principle is inspired by the work
of Mittnik and Paolella (2000) and Robben et al. (2022), and is motivated by
the fact that earlier years’ mortality dynamics are less representative
compared to the dynamics observed in more recent years. This is in particular
true when a lengthy calibration period is used as in our case study. The
estimated mean parameter $\hat{\boldsymbol{c}}$ and covariance matrix
$\hat{\boldsymbol{\Sigma}}_{w}$ equal:

$\displaystyle\hat{\boldsymbol{c}}=\begin{bmatrix}[r]0.01170\\\ -0.06590\\\
0\\\ 0\end{bmatrix},\hskip
14.22636pt\hat{\boldsymbol{\Sigma}}_{w}=\begin{bmatrix}[r]0.00495&0.00012&0.00052&-0.00120\\\
0.00012&0.01491&0.00281&-0.00320\\\ 0.00052&0.00281&0.00539&-0.00321\\\
-0.00120&-0.00320&-0.00321&0.01564\end{bmatrix}.$

We generate $10\ 000$ trajectories for the vector of period effects
$\boldsymbol{\mathcal{K}}_{t}$ for $t\in\\{1851,...,2021\\}$ by taking random
draws from a four-dimensional Gaussian distribution with the estimated
$\hat{\boldsymbol{c}}$ and $\hat{\boldsymbol{\Sigma}}_{w}$ as parameters for
the mean and covariance matrix, respectively. Figure 12 shows the $0.5\%$,
median and $99.5\%$ quantile, based on these $10\ 000$ trajectories for
$K_{t}^{(1)}$, $K_{t}^{(2)}$, $\kappa_{t}^{(1)}$ and $\kappa_{t}^{(2)}$.

Figure 12: The $0.5\%$, median and $99.5\%$ quantile, based on $10\ 000$
trajectories for the estimated time series models of the four-dimensional
vector of period effects $\hat{\boldsymbol{\mathcal{K}}}_{t}$ (blue fan
chart). The blue dots represent the calibrated period effects. The green line
represents one trajectory for each period effect. Projection period 2022-2080.

####  6.4.2 Projecting the regime-switching model to generate age-specific
mortality shocks

We now focus on generating the mortality shocks for the age groups
$\mathcal{X}^{(1)}$, i.e. ages 20-59, and $\mathcal{X}^{(2)}$, i.e. ages
60-85. To better align the projections with real-world settings, we suggest
two minor modifications to the projection strategy outlined in Section 5.2.
These modifications are explained in detail in Appendix E.2 and ensure that
the generated shocks in a high volatility period do not have a long-lasting
impact on future mortality rate projections. We then generate 10 000
trajectories for both the Markov chains and the regime-switching models,
covering the period from 2022 to 2080. Figure 13 displays one such trajectory.
First, related to the COVID-19 pandemic, we assume that we leave the high
volatility regime in the year 2023 for the age group 60-85.262626Related to
this, we also make sure that the mortality shock related to COVID-19 in the
years 2020-2022 has been offset in 2023, see Appendix E.2. We further observe
two future spells in the high volatility regime of three and two years
respectively in which mortality shocks occur.

Figure 13: In the left panel, we visualize one trajectory of the Markov chain
$\rho_{t}^{(2)}$ for the age group 60-85 over the projection period 2022-2080.
In the right panel, we show the corresponding generated trajectory for the
regime-switching model with age-specific effect
$\hat{\mathfrak{B}}_{x}Y_{t}^{(2)}$ averaged over the ages in the interval
60-85. The underlying Markov chain occupies the states pictured in the left
panel. The light blue vertical bars refer to periods in the high volatility
regime.

####  6.4.3 Mortality rate projections

Following Section 5.2, we forecast the force of mortality for an $x$-year old
individual in year $t\in\mathcal{T}^{\text{{\tiny
pred}}}:=\\{2022,...,2080\\}$ recursively as

$\hat{\mu}_{x,t,\iota}^{({\scriptscriptstyle\text{NL}})}=\hat{\mu}_{x,t-1,\iota}^{({\scriptscriptstyle\text{NL}})}\cdot\exp\left(\hat{A}_{x}+\displaystyle\sum_{i=1}^{2}\hat{B}_{x}^{(i)}\hat{K}_{t,\iota}^{(i)}+\displaystyle\sum_{j=1}^{2}\hat{\beta}_{x}^{(j)}\hat{\kappa}_{t,\iota}^{(j)}+\hat{\mathfrak{B}}_{x}^{(a)}\hat{Y}_{t,\iota}^{(a)}\right),$

for age group $a\in\\{[20,59],[60,85]\\}$, for $t>2021$ and referring to the
$\iota$-th generated trajectory. As starting value of the recursion in the
year $2021$ we use the observed crude central death rate, i.e.
$\hat{\mu}_{x,2021,\iota}^{({\scriptscriptstyle\text{NL}})}=d_{x,2021}^{({\scriptscriptstyle\text{NL}})}/E_{x,2021}^{({\scriptscriptstyle\text{NL}})}$
for all $\iota$. Using Equation (2.2), we can calculate the corresponding
mortality rates.

Figure 14 illustrates three trajectories of the mortality rates at ages 35
(left panel) and 80 (right panel) over the projection period 2022-2080. The
generated mortality decline is occasionally disrupted by some artificially
generated mortality shocks. Multiple mortality shocks, sometimes with very
high severity, occur in the left panel of Figure 14. This is mainly due to the
integration of the two World Wars in the calibration procedure of the regime-
switching model which mainly affects younger age groups. Section 6.5 further
elaborates on this.

Figure 14: Three generated trajectories for the mortality rates at ages 35 and
80. Calibration period 1850-2021 and projection period 2022-2080. We zoom in
on the years starting from 1950 to better visualize the projections.

Figure 15 shows the projected mortality rates for individuals aged 35, 50, 65,
and 80. The plot displays the $90\%$, $95\%$, and $99\%$ prediction intervals
based on $10\ 000$ generated trajectories. The blue dots indicate the observed
mortality rates and the blue line represents the fit of the baseline mortality
improvement model. The green line shows the in-sample fit of the mortality
model constructed by the Dutch Actuarial Association in 2020, i.e. the AG2020
model, a Li-Lee multi-population mortality model calibrated over the 1970-2019
period (Royal Dutch Actuarial Association, 2020). The green fan charts display
the $90\%$, $95\%$, and $99\%$ prediction intervals based on the simulations
resulting from the AG2020 model. Our methodology results in wider fan charts
than the AG2020 model due to the allowance for future mortality shocks. The
median quantile of the simulations resulting from our proposed regime-
switching mortality improvement model and the one resulting from the AG2020
model approximately coincide. As a result of the modifications made to the
projection strategy described in Appendix E.2, the bottom right panel of
Figure 15 demonstrates that the COVID-19 shock has been offset by the year
2022 for individuals aged 80.

Figure 15: The blue fan charts visualize the $90\%$, $95\%$, and $99\%$
prediction intervals based on $10\ 000$ generated trajectories of the
mortality improvement model with shock regime. We use calibration period
1850-2021, age range 20-85, and projection period 2022-2080. The blue shaded
dots are the observed mortality rates and the blue line shows the fit of the
baseline mortality improvement model. Additionally, the green line shows the
in-sample fit of the AG2020 model where we calibrated the Li-Lee model on data
from 1970-2019. The green fan charts show the $90\%$, $95\%$, and $99\%$
prediction intervals based on the simulations resulting from this AG2020
model. The median quantiles of the simulations for both methods are shown in
dark blue and dark green, respectively.

###  6.5 Sensitivity analysis

To reflect differences in risk appetite, we perform a sensitivity analysis
with our proposed toolbox. The inclusion of World War I and World War II has a
significant impact on the calibration of the regime-switching model. Figure 9
shows that the largest outliers in the empirical distribution of the residuals
$Z_{x,t}$ occur during the periods of the two world wars. This consequently
leads to a large value for the calibrated volatility parameter
$\hat{\sigma}_{H}$ in the high volatility regime, see Table 1. The wide
confidence intervals observed in the projected mortality rates at ages 35 and
50 in Figure 15 confirm this issue. Furthermore, insurance contracts nowadays
often include a war clause that excludes coverage during wars. Therefore, the
risk modeler may want to exclude these war events when calibrating the regime-
switching model.

We explore a specific scenario that calibrates the improvement model on data
from 1875, eliminating war events such as the Crimean and Franco-German wars
that happened in the period 1850-1875, see Figure 2. Furthermore, we exclude
World War I and II in the calibration of the shock regime by introducing
weights in its log-likelihood function. Equation (4.7) becomes:

$\displaystyle l(\Theta)$ $\displaystyle=\sum_{t\in\mathcal{T}}\nu_{t}\log
f\left(\boldsymbol{z}_{t}\mid\boldsymbol{z}_{t-1},...,\boldsymbol{z}_{t_{\min}},\Theta\right),$
(6.8)

and we choose $\nu_{t}=0$ for
$t\in\\{1914\text{-}1919,1940\text{-}1946\\}$.272727Note that we also exclude
the years 1919 and 1946 since we work with a mortality improvement model and
the years 1919 and 1946 are as such impacted by the world wars as well. The
calibration of the regime-switching model then occurs in the same way as
explained in Section 4.2 with details in Appendix E.

Table 2 and Figure 16 display the calibrated parameters and age-specific
effects $\hat{\mathfrak{B}}_{x}$ in the regime-switching model for the age
ranges 20-59 and 60-85. Excluding war events leads to a significant decrease
in $\sigma_{H}$, the volatility parameter in the high volatility regime. For
the 20-59 age range, $\sigma_{H}$ drops from 2.41 to 0.62, and for the 60-85
age range, it drops from 0.90 to 0.39. The drop in $\sigma_{H}$ for the 20-59
age range is due to the exclusion of the war victims who died at a young age.
Excluding all war events also reduces the probability of transitioning from
the low to the high volatility regime, from 0.047 to 0.019 for the 20-59 age
range, and from 0.67 to 0.052 for the 60-85 age range. The values of the
calibrated volatility and slope parameters in the low volatility regime
display only minor changes, which confirms the stability of the calibration
strategy.

Parameter 20-59 60-85 $p_{1,2}$ 0.01873 0.05239 $p_{2,1}$ 0.37218 0.78693 $\sigma_{e_{1}}$ 0.12612 0.05394 $\text{slope}_{1}$ -0.00171 0.00050 $\sigma_{e_{2}}$ 0.18224 0.03919 $\text{slope}_{2}$ -0.00344 -0.00050 $\mu_{H}$ 0.02684 -0.01665 $\sigma_{H}$ 0.61850 0.39284 |  x
---|---
Table 2: Calibrated parameters in the regime-switching model for the two age groups. We exclude the world wars during the calibration. Calibration period 1875-2021, male data. |  Figure 16: Calibrated age-specific effect $\mathfrak{B}_{x}$ in the regime-switching model for the two age groups. We exclude the world wars during the calibration. Calibration period 1875-2021, male data.

Figure 17 displays the fan charts that represent the $90\%$, $95\%$, and
$99\%$ prediction intervals based on the projected mortality rates spanning
the years 2022-2080 for individuals aged 35, 50, 65, and 80. Notably, in
contrast to Figure 15, the variability in the generated mortality rates
significantly reduces, but still mostly exceeds the variability from the
AG2020 model.

Figure 17: The blue fan charts visualize the $90\%$, $95\%$, and $99\%$
prediction intervals based on $10\ 000$ generated trajectories constructed
with the mortality improvement model with shock regime. We use calibration
period 1875-2021, age range 20-85, and projection period 2022-2080. We further
exclude the two world wars during the calibration of the regime-switching
model. The blue shaded dots are the observed mortality rates and the blue line
shows the fit of the baseline mortality improvement model. Additionally, the
green line shows the in-sample fit of the AG2020 model where we calibrated the
Li-Lee model on data from 1970-2019. The green fan charts show the $90\%$,
$95\%$, and $99\%$ prediction intervals based on the simulations resulting
from this AG2020 model. The median quantiles of the simulations for both
methods are shown in dark blue and dark green, respectively.

###  6.6 Solvency capital requirement

We focus on the calculation of the solvency capital requirement (SCR) for the
life underwriting risk module in the Solvency II Directive (2009/138/EC), and
examine the submodules of mortality, longevity, and catastrophe risk. We
explore two approaches: the run-off VaR, as discussed by Gylys and Šiaulys
(2019), Richards et al. (2014) and Levantesi and Pizzorusso (2019), and the
standard model outlined by European Commission (2010). The standard model,
unlike the VaR approach, does not require stochastic modeling of mortality
rates. The run-off VaR relies on mortality rates generated from a stochastic
model for the future years of the insurance coverage.282828Alternatively, one
could consider the one-year VaR which uses stochastically simulated mortality
rates only for the first projection year. One approach is to recalibrate the
mortality model based on the outcomes of the first-year simulations, and
relies on the projected best-estimate mortality rates of these recalibrated
models for the subsequent years (Richards et al., 2014). See, e.g., Gylys and
Šiaulys (2019) for more details and a comparison between the one-year VaR and
the run-off VaR. We hereby use the scenarios generated from our proposed
mortality improvement model with shock regime. Appendix F provides more
details about the methodologies and their implementation. We put focus on the
valuation of an immediate life annuity in Section 6.6.1 that is solely exposed
to longevity risk and on the valuation of a term life insurance product in
Section 6.6.2, which is exposed to mortality and catastrophe risk.

####  6.6.1 SCR for an immediate life annuity

We consider the valuation of an immediate whole life annuity for an individual
aged $x$ in the year 2021, with an annual payout at the end of each year of
€10 000 until the insured dies. An immediate life annuity is exposed to
longevity risk since the insurer needs to pay out more annuities if the
policyholder lives longer than expected. The standard model obtains the SCR
associated with longevity risk as the change in Best Estimate Liabilities
(BEL) under a longevity shock of $20\%$. In this shock scenario, all future
mortality rates are multiplied by a factor of 0.8. The BEL in the year 2021
reflects the present value of future cash flows from the insurer to the
policyholder, discounted to the year 2021, and is determined based on best-
estimate projections of mortality rates and survival probabilities. To
determine the run-off VaR, we use our proposed mortality improvement model and
generate 10 000 simulations of the expected present value of the insurer’s
liabilities in the year 2021. Subsequently, we calculate the $99.5\%$ quantile
of these simulated liabilities and subtract the BEL in the year 2021. Appendix
F.1 further details the calculations of the SCR using both the run-off VaR and
the standard model approach.

The left panel of Figure 18 presents the SCR for an immediate life annuity
starting in 2021 for individuals aged between 55 and 90. The SCR obtained with
the standard model first slightly increases and then decreases faster with
age. Moreover, the discrepancy between the SCR derived through the VaR
approach and the SCR obtained with the standard model is more pronounced for
older ages. This phenomenon can be attributed to the fact that the impact of a
20$\%$ decrease in mortality rates on the Best Estimate Liabilities (BEL) is
larger for older age mortality rates compared to younger age mortality rates,
see Equation (F.5). These findings are consistent with the results obtained in
Börger (2010). Therefore, a longevity portfolio predominantly comprising older
individuals may be overly cautious regarding the amount of solvency capital
when implementing the standard model. Instead of a fixed, age-independent
longevity shock on all mortality rates, an age-dependent longevity shock would
bring the results obtained with the standard model and the VaR approach closer
to each other. The findings from the longevity stress test in EIOPA (2018),
who derive age-dependent longevity shocks, support these results. However,
EIOPA recommends retaining the age-independent longevity shock in the standard
model, citing concerns about added complexity and implementation costs.

Figure 18: In the left panel, the SCR for an immediate life annuity contract
paying €10 000 annually and starting in the year 2021 for a person aged 55-90.
In the right panel, the SCR for a term life insurance contract with a death
benefit of €100 000 starting in the year 2021 for a person aged 20-60, ending
in the year the policyholder reaches the age 65. The dark blue line shows the
SCR calculated with the run-off VaR approach and the light blue line is the
one calculated with the standard model.

####  6.6.2 Term life insurance

We examine a term life insurance policy with a terminal age of 65 that is
issued to an individual aged $x$ in the year 2021. The insured death benefit
equals €150 000 and is payable at the end of the year of death. A term life
insurance is exposed to mortality and catastrophe risk as the probability of
paying out the death benefit before reaching the term of the contract
increases when the mortality rates increase. The standard model obtains the
SCR associated with mortality risk as the change in the BEL under an upward
mortality shock of $15\%$ (European Commission, 2010). The SCR associated with
catastrophe risk, under the standard model, is obtained as the change in BEL
under an absolute increase in the mortality rates of the year 2021 by 0.0015.
These two separate SCRs are aggregated according to the guidelines in European
Commission (2010). Next, similarly to Section 6.6.1, we obtain the SCR with
the run-off VaR approach by generating $10\ 000$ simulations for future
mortality rates including the occurrence of shocks. Appendix F.2 provides the
relevant details.

The right panel of Figure 18 shows the SCR for a term life insurance issued in
the year 2021 to an individual aged between 20 and 60. For younger ages, i.e.
for term insurance policies with a longer term, the SCR obtained with the
standard model is lower than the SCR using our proposed mortality model.
Conversely, for ages above 45 years, we observe the opposite. This stems from
two factors: the increased impact of a $15\%$ mortality shock on the BEL for
older age mortality rates in the standard model and the greater uncertainty in
mortality rate predictions for younger ages compared to older ages in the VaR
approach, see Figure 17. In line with our findings in Section 6.6.1, we
emphasize the need for incorporating an age-dependent mortality and
catastrophe shock in the standard model.

##  7 Discussion

In this paper, we have developed a framework for a stochastic multi-population
mortality projection model that incorporates the occurrence and severity of
age-specific mortality shocks into future mortality projections using a
regime-switching model. We showed that the width of the mortality prediction
intervals largely depends on which mortality shocks are considered in the
calibration procedure of the regime-switching model.

From a technical point of view, we discussed the importance of mortality
improvement models as an alternative to the classic mortality models used in
most academic literature. Overall, the proposed model can assist in modeling
future mortality rates, which is crucial for actuaries, life insurers, and
risk managers in light of the need to adhere to the capital requirements
imposed by supervisory authorities. In terms of future research directions, a
possible avenue is to extend the current two-state regime model to a three-
state regime model. In this setting, regime 1 will represent slow mortality
improvements, regime 2 fast mortality improvements, and regime 3 will capture
mortality shocks. This extension incorporates a possible switch to slower
mortality improvements, a behavior that has been observed in the last decade
(Djeundje et al., 2022). We demonstrated the calibration of the regime-
switching model under the assumption of normally distributed shocks. Future
research can focus on investigating alternative distributional assumptions,
their impact on the generated scenarios, and the development of a suitable
calibration strategy for the regime-switching model.

From an economic point of view, our paper contributes to the literature on
scenario thinking in risk management, which considers different hypothetical
future scenarios and their corresponding implications. Hereto, we perform a
quantitative impact assessment that investigates the impact of different
scenarios related to the inclusion or exclusion of mortality shocks into
future mortality rate projections. Using these different scenarios,
(re)insurers can assess the corresponding impact on the calculation of, e.g.,
solvency capital requirements, and tailor the proposed framework to their risk
management objectives. Scenario thinking is gaining attention in the actuarial
literature. Gielen and de Waegenaere (2014) consider specific scenarios
involving a 100$\%$, 50$\%$, or $20\%$ reduction in the one-year mortality
numbers caused by cancer. They then evaluate the corresponding impact on
capital requirements in comparison to traditional longevity models such as the
Lee-Carter model. Hanika (2023) studies the combined effects of a financial
and a mortality shock on a life insurer’s risk during a COVID-19 stress test.
He emphasizes the relevance of scenario thinking by evaluating practical risk
mitigation strategies for mitigating the short-term impact of pandemic events,
offering insights for future pandemic scenarios. Bongiorno et al. (2022)
exemplify how pension schemes can incorporate climate scenarios into their
considerations, including aspects like funding strategies and long-term
strategic asset allocation.

## References

  * Antonio et al. (2017) Katrien Antonio, Sander Devriendt, Wouter de Boer, Robert de Vries, Anja De Waegenaere, Hok-Kwan Kan, Egbert Kromme, Wilbert Ouburg, Tim Schulteis, Erica Slagter, et al.  Producing the Dutch and Belgian mortality projections: a stochastic multi-population standard.  _European Actuarial Journal_ , 7(2):297–336, 2017\.  doi: 10.1007/s13385-017-0159-x.
  * Antonio et al. (2020) Katrien Antonio, Sander Devriendt, and Jens Robben.  The IA$\mid$BE 2020 mortality projection for the Belgian population, 2020.  Available online: https://katrienantonio.github.io/papers/report_iabe_2020.pdf.
  * Bongiorno et al. (2022) Luca Bongiorno, Andrew Claringbold, Lisa Eichler, Claire Jones, Bert Kramer, Louise Pryor, and Nick Spencer.  Climate scenario analysis for pension schemes: a UK case study.  _British Actuarial Journal_ , 27:e6, 2022.  doi: 10.1017/S1357321721000106.
  * Börger (2010) Matthias Börger.  Deterministic shock vs. stochastic Value–at–Risk—an analysis of the Solvency II standard model approach to longevity risk.  _Blätter der DGVFM_ , 31:225–259, 2010.  doi: 10.1007/s11857-010-0125-z.
  * Brest et al. (2006) Janez Brest, Sao Greiner, Borko Boskovic, Marjan Mernik, and Viljem Zumer.  Self-adapting control parameters in differential evolution: A comparative study on numerical benchmark problems.  _IEEE Transactions on Evolutionary Computation_ , 10(6):646–657, 2006.  doi: 10.1109/TEVC.2006.872133.
  * Brouhns et al. (2002) Natacha Brouhns, Michel Denuit, and Jeroen K Vermunt.  A Poisson log-bilinear regression approach to the construction of projected lifetables.  _Insurance: Mathematics and Economics_ , 31(3):373–393, 2002.  doi: 10.1016/S0167-6687(02)00185-3.
  * Cairns et al. (2009) Andrew JG Cairns, David Blake, Kevin Dowd, Guy D Coughlan, David Epstein, Alen Ong, and Igor Balevich.  A quantitative comparison of stochastic mortality models using data from England and Wales and the United States.  _North American Actuarial Journal_ , 13(1):1–35, 2009.  doi: 10.1080/10920277.2009.10597538.
  * Chen and Cox (2009) Hua Chen and Samuel H. Cox.  Modeling mortality with jumps: Applications to mortality securitization.  _The Journal of Risk and Insurance_ , 76(3):727–751, 2009.  doi: 10.1111/j.1539-6975.2009.01313.x.
  * Chen and Cummins (2010) Hua Chen and J. David Cummins.  Longevity bond premiums: The extreme value approach and risk cubic pricing.  _Insurance: Mathematics and Economics_ , 46(1):150–161, 2010.  doi: 10.1016/j.insmatheco.2009.09.007.
  * Cleveland and Tiao (1976) William P Cleveland and George C Tiao.  Decomposition of seasonal time series: A model for the census X-11 program.  _Journal of the American Statistical Association_ , 71(355):581–587, 1976.  doi: 10.1080/01621459.1976.10481532.
  * Continuous Mortality Investigation Bureau (2002) Continuous Mortality Investigation Bureau.  An interim basis for adjusting the “92” series mortality projections for cohort effects, 2002.  Available online: https://www.actuaries.org.uk/system/files/documents/pdf/cmiwp1.pdf.
  * Cox et al. (2006) Samuel H. Cox, Yijia Lin, and Shaun Wang.  Multivariate exponential tilting and pricing implications for mortality securitization.  _Journal of Risk and Insurance_ , 73(4):719–736, 2006.  doi: 10.1111/j.1539-6975.2006.00196.x.
  * Currie (2006) Iain D Currie.  Smoothing and forecasting mortality rates with P-splines, 2006.  Talk given at the Institute of Actuaries. Available online: https://www.actuaries.org.uk/system/files/documents/pdf/currie.pdf.
  * Djeundje et al. (2022) Viani B Djeundje, Steven Haberman, Madhavi Bajekal, and Joseph Lu.  The slowdown in mortality improvement rates 2011–2017: a multi-country analysis.  _European Actuarial Journal_ , 12(2):839–878, 2022.  doi: 10.1007/s13385-022-00318-0.
  * Dowd et al. (2020) Kevin Dowd, Andrew JG Cairns, and David Blake.  CBDX: a workhorse mortality model from the Cairns–Blake–Dowd family.  _Annals of Actuarial Science_ , 14(2):445–460, 2020.  doi: 10.1017/S1748499520000159.
  * EIOPA (2018) EIOPA.  EIOPA’s second set of advice to the European Commission on specific items in the Solvency II delegated regulation, 2018.  Available online: https://www.eiopa.europa.eu/system/files/2019-04/eiopa-18-075-eiopa_second_set_of_advice_on_sii_dr_review.pdf.
  * European Commission (2010) European Commission.  QIS5 technical specifications, 2010.  Available online: https://register.eiopa.europa.eu/Publications/QIS/QIS5TS20100706.doc.
  * Gielen and de Waegenaere (2014) Jeroen Gielen and Anja de Waegenaere.  What if there was a cure for cancer? Improving the quantification of longevity risk using scenarios.  _The European Actuary_ , 5(1):10–11, 2014.
  * Goes et al. (2023) Julius Goes, Karim Barigou, and Anne Leucht.  Bayesian mortality modelling with pandemics: a vanishing jump approach, 2023.  Available on arXiv: doi: 10.48550/arXiv.2311.04920.
  * Golub et al. (1979) Gene H Golub, Michael Heath, and Grace Wahba.  Generalized cross–validation as a method for choosing a good ridge parameter.  _Technometrics_ , 21(2):215–223, 1979.  doi: 10.2307/1268518.
  * Gungah and Narsoo (2021) Gavakshi Gungah and Jason Narsoo.  A novel EVT-modified Lee-Carter model for mortality forecasting: An application to extreme mortality events.  _Journal of Statistics and Management Systems_ , 25(1):211–243, 2021.  doi: 10.1080/09720510.2021.1878629.
  * Gylys and Šiaulys (2019) Rokas Gylys and Jonas Šiaulys.  Revisiting calibration of the Solvency II standard formula for mortality risk: Does the standard stress scenario provide an adequate approximation of Value-at-Risk?  _Risks_ , 7(2), 2019.  doi: 10.3390/risks7020058.
  * Haberman and Renshaw (2011) Steven Haberman and Arthur Renshaw.  A comparative study of parametric mortality projection models.  _Insurance: Mathematics and Economics_ , 48(1):35–55, 2011.  doi: 10.1016/j.insmatheco.2010.09.003.
  * Haberman et al. (2014) Steven Haberman, Vladimir Kaishev, Pietro Millossovich, Andrés Villegas, Steven Baxter, Andrew Gaches, Sveinn Gunnlaugsson, and Mario Sison.  Longevity basis risk: A methodology for assessing basis risk.  Technical report, Institute and Faculty of Actuaries, 2014.  Available online: https://www.actuaries.org.uk/system/files/documents/pdf/ifoa-llma-longevity-basis-risk-report.pdf.
  * Hainaut (2012) Donatien Hainaut.  Multidimensional Lee-–Carter model with switching mortality processes.  _Insurance: Mathematics and Economics_ , 50(2):236–246, 2012.  doi: 10.1016/j.insmatheco.2011.11.003.
  * Hamilton (1989) James D Hamilton.  A new approach to the economic analysis of nonstationary time series and the business cycle.  _Econometrica_ , 57(2):357–384, 1989.  doi: 10.2307/1912559.
  * Hanika (2023) Moritz Hanika.  A COVID-19 stress test for life insurance: insights into the effectiveness of different risk mitigation strategies.  _European Actuarial Journal_ , pages 1–26, 2023.  doi: 10.1007/s13385-023-00371-3.
  * Hardy (2001) Mary R Hardy.  A regime–switching model of long–term stock returns.  _North American Actuarial Journal_ , 5(2):41–53, 2001.  doi: 10.1080/10920277.2001.10595984.
  * Hunt and Blake (2020) Andrew Hunt and David Blake.  Identifiability in age/period mortality models.  _Annals of Actuarial Science_ , 14(2):461––499, 2020.  doi: 10.1017/S1748499520000111.
  * Hunt and Villegas (2021) Andrew Hunt and Andrés M Villegas.  Mortality improvement rates: Modeling, parameter uncertainty, and robustness.  _North American Actuarial Journal_ , pages 1–27, 2021.  doi: 10.1080/10920277.2021.2006068.
  * Kannisto (1994) Väinö Kannisto.  _Development of oldest-old mortality, 1950-1990: Evidence from 28 developed countries_ , volume 1.  Odense University Press, 1994.
  * Lee (2003) Ronald Lee.  The demographic transition: three centuries of fundamental change.  _Journal of economic perspectives_ , 17(4):167–190, 2003.  doi: 10.1257/089533003772034943.
  * Lee and Carter (1992) Ronald D Lee and Lawrence R Carter.  Modeling and forecasting U.S. mortality.  _Journal of the American Statistical Association_ , 87(419):659–671, 1992.  doi: 10.2307/2290201.
  * Levantesi and Pizzorusso (2019) Susanna Levantesi and Virginia Pizzorusso.  Application of machine learning to mortality modeling and forecasting.  _Risks_ , 7(1):26, 2019.  doi: 10.3390/risks7010026.
  * Li and Lee (2005) Nan Li and Ronal Lee.  Coherent mortality forecasts for a group of populations: an extension of the Lee–Carter method.  _Demography_ , 42(3):575–594, 2005.  doi: 10.1353/dem.2005.0021.
  * Li and Chan (2005) Siu-Hang Li and Wai-Sum Chan.  Outlier analysis and mortality forecasting: The United Kingdom and Scandinavian countries.  _Scandinavian Actuarial Journal_ , 2005(3):187–211, 2005.  doi: 10.1080/03461230510006973.
  * Liu and Li (2015) Yanxin Liu and Johnny Siu-Hang Li.  The age pattern of transitory mortality jumps and its impact on the pricing of catastrophic mortality bonds.  _Insurance: Mathematics and Economics_ , 64:135–150, 2015\.  doi: 10.1016/j.insmatheco.2015.05.005.
  * Milidonis et al. (2011) Andreas Milidonis, Yijia Lin, and Samuel H Cox.  Mortality regimes and pricing.  _North American Actuarial Journal_ , 15(2):266–289, 2011.  doi: 10.1080/10920277.2011.10597621.
  * Mitchell et al. (2013) Daniel Mitchell, Patrick Brockett, Rafael Mendoza-Arriaga, and Kumar Muthuraman.  Modeling and forecasting mortality rates.  _Insurance: Mathematics and Economics_ , 52(2):275–285, 2013.  doi: 10.1016/j.insmatheco.2013.01.002.
  * Mittnik and Paolella (2000) Stefan Mittnik and Marc S Paolella.  Conditional density and Value–at–Risk prediction of Asian currency exchange rates.  _Journal of Forecasting_ , 19(4):313–333, 2000\.  doi: 10.1002/1099-131X(200007)19:4¡313::AID-FOR776¿3.0.CO;2-E.
  * Moritz and Bartz-Beielstein (2017) Steffen Moritz and Thomas Bartz-Beielstein.  imputeTS: time series missing value imputation in R.  _The R Journal_ , 9(1):207, 2017.  doi: 10.32614/RJ-2017-009.
  * Pitacco et al. (2009) Ermanno Pitacco, Michel Denuit, Steven Haberman, and Annamaria Olivieri.  _Modeling Longevity Dynamics for Pensions and Annuity Business_.  Oxford University Press, London, 2009.
  * Plat (2009) Richard Plat.  On stochastic mortality modeling.  _Insurance: Mathematics and Economics_ , 45(3):393–404, 2009.  doi: 10.1016/j.insmatheco.2009.08.006.
  * Price et al. (2006) Kenneth Price, Rainer M Storn, and Jouni A Lampinen.  _Differential evolution: a practical approach to global optimization_.  Springer Science & Business Media, 2006.
  * Renshaw and Haberman (2003) Arthur E Renshaw and Steven Haberman.  Lee–Carter mortality forecasting with age-specific enhancement.  _Insurance: Mathematics and Economics_ , 33(2):255–272, 2003.  doi: 10.1016/S0167-6687(03)00138-0.
  * Renshaw and Haberman (2006) Arthur E Renshaw and Steven Haberman.  A cohort–based extension to the Lee–Carter model for mortality reduction factors.  _Insurance: Mathematics and Economics_ , 38(3):556–570, 2006.  doi: 10.1016/j.insmatheco.2005.12.001.
  * Richards et al. (2014) Stephen J Richards, Ian D Currie, and Gavin P Ritchie.  A Value–at–Risk framework for longevity trend risk.  _British Actuarial Journal_ , 19(1):116–139, 2014\.  doi: 10.1017/S1357321712000451.
  * Robben et al. (2022) Jens Robben, Katrien Antonio, and Sander Devriendt.  Assessing the impact of the COVID-19 shock on a stochastic multi-population mortality model.  _Risks_ , 10(2), 2022.  doi: 10.3390/risks10020026.
  * Royal Dutch Actuarial Association (2020) Royal Dutch Actuarial Association.  Projections life table AG2020, 2020.  Available online: https://www.actuarieelgenootschap.nl/download/projections-life-table-ag2020-pdf.htm.
  * Royal Dutch Actuarial Association (2022) Royal Dutch Actuarial Association.  Projections life table AG2022, 2022.  Available online: https://www.actuarieelgenootschap.nl/download/projections-life-table-ag2022.htm.
  * Toole (2007) Jim Toole.  Potential impact of pandemic influenza on the US life insurance industry.  _Society of Actuaries_ , 2007.
  * Zhou et al. (2013) Rui Zhou, Johnny Siu-Hang Li, and Ken Seng Tan.  Pricing standardized mortality securitizations: A two-population model with transitory jump effects.  _Journal of Risk and Insurance_ , 80(3):733–774, 2013.  doi: 10.1111/j.1539-6975.2013.12015.x.

##  Appendix A Identifiability constraints in a two-factor mortality
improvement model

##### Model specifications.

Renshaw and Haberman (2003) introduced a two-factor version of the Lee-Carter
model. We adopt a similar specification for the mortality improvement rates,
as follows:

$\displaystyle\log\mu_{x,t}-\log\mu_{x,t-1}=A_{x}+B_{x}^{(1)}K_{t}^{(1)}+B_{x}^{(2)}K_{t}^{(2)},$
(A.1)

where $\mu_{x,t}$ denotes the force of mortality at age $x\in\mathcal{X}$ and
year $t\in\mathcal{T}$. We denote the age-specific effects as
$A_{x},\>B_{x}^{(1)}$ and $B_{x}^{(2)}$ and the period-specific effects as
$K_{t}^{(1)}$ and $K_{t}^{(2)}$. We introduce vector notations for these age
and period effects as follows:

$\boldsymbol{B}^{(i)}=\left(B_{x}^{(i)}\right)_{x\in\mathcal{X}}\in\mathbb{R}^{|\mathcal{X}|\times
1},\hskip
14.22636pt\boldsymbol{K}^{(i)}=\left(K_{t}^{(i)}\right)_{t\in\mathcal{T}}\in\mathbb{R}^{|\mathcal{T}|\times
1}\hskip 21.33955pt\text{for}\>i\in\\{1,2\\},$

where $|\mathcal{X}|$ and $|\mathcal{T}|$ denote the number of ages and years
considered in the age range $\mathcal{X}$ and calibration period
$\mathcal{T}$, respectively.

##### Introducing the parameter constraints.

To make the mortality improvement model in Equation (A.1) fully identifiable,
we need to impose parameter constraints. Straightforward extensions of the
identifiability constraints in the Lee-Carter model are:

$||\boldsymbol{B}^{(i)}||_{2}^{2}=\displaystyle\sum_{x\in\mathcal{X}}\left(B_{x}^{(i)}\right)^{2}=1,\hskip
14.22636pt\sum_{t\in\mathcal{T}}K_{t}^{(i)}=0,\hskip
28.45274pt\text{for}\>i\in\\{1,2\\},$ (A.2)

where $||\cdot||_{2}$ denotes the Euclidean norm. We use constraints on the
squared age effects $B_{x}^{(i)}$ for $i\in\\{1,2\\}$ to enhance numerical
stability. However, these constraints alone are not sufficient to fully
identify the model in Equation (A.1). For instance, some transformations of
the age and period effects can still lead to the same model fit:

$\displaystyle A_{x}$ $\displaystyle\mapsto A_{x},$ (A.3) $\displaystyle
B_{x}^{(1)}$ $\displaystyle\mapsto B_{x}^{(1)}+B_{x}^{(2)},\hskip
14.22636ptB_{x}^{(2)}\mapsto B_{x}^{(2)},$ $\displaystyle K_{t}^{(1)}$
$\displaystyle\mapsto K_{t}^{(1)},\hskip 45.52458ptK_{t}^{(2)}\mapsto
K_{t}^{(2)}-K_{t}^{(1)}.$

We refer to such a transformation as an invariant transformation. Hunt and
Blake (2020) note that for a Lee-Carter model with $N$ period effects,
$N\cdot(N+1)$ parameter constraints are needed. We hence need an additional
two constraints to fully identify the mortality improvement model. A
convenient and useful set of constraints are the orthogonality constraints on
the age and period effects (Hunt and Blake, 2020), i.e. we impose:

$\displaystyle\langle\boldsymbol{B}^{(1)},\boldsymbol{B}^{(2)}\rangle=\displaystyle\sum_{x\in\mathcal{X}}B_{x}^{(1)}B_{x}^{(2)}=0,\hskip
14.22636pt\langle\boldsymbol{K}^{(1)},\boldsymbol{K}^{(2)}\rangle=\displaystyle\sum_{t\in\mathcal{T}}K_{t}^{(1)}K_{t}^{(2)}=0,$
(A.4)

where $\langle\cdot,\cdot\rangle$ refers to the Euclidean inner product. This
makes transformations as in Equation (A.3) impossible. In the remainder of
this appendix, we demonstrate how to convert a set of parameters
$(A_{x},B_{x}^{(1)},B_{x}^{(2)},K_{t}^{(1)},K_{t}^{(2)})$ into a set of
parameters
$(\tilde{A}_{x},\tilde{B}_{x}^{(1)},\tilde{B}_{x}^{(2)},\tilde{K}_{t}^{(1)},\tilde{K}_{t}^{(2)})$
that satisfy the identifiability constraints defined in Equations (A.2) and
(A.4), and result in equivalent model fits, see Equation (A.1).

##### Linear transformations.

Hunt and Blake (2020) show that the only possible invariant transformations
are of the form:

$\displaystyle\tilde{A}_{x}$
$\displaystyle:=A_{x}-\boldsymbol{B}_{x}^{T}\boldsymbol{A}\boldsymbol{D},$
(A.5) $\displaystyle\tilde{\boldsymbol{B}}_{x}^{T}$
$\displaystyle=\left[\tilde{B}_{x}^{(1)},\tilde{B}_{x}^{(2)}\right]^{T}:=\boldsymbol{B}_{x}^{T}\boldsymbol{A}$
$\displaystyle\tilde{\boldsymbol{K}}_{t}$
$\displaystyle=\left[\tilde{K}_{t}^{(1)},\tilde{K}_{t}^{(2)}\right]:=\boldsymbol{A}^{-1}\boldsymbol{K}_{t}+\boldsymbol{D},$

where $\boldsymbol{B}_{x}=(B_{x}^{(1)},B_{x}^{(2)})\in\mathbb{R}^{2\times 1}$,
T denotes the transpose operator, $\boldsymbol{A}\in\mathbb{R}^{2\times 2}$ is
invertible and $\boldsymbol{D}\in\mathbb{R}^{2\times 1}$. Note that the above
transformations contain six free parameters, i.e. four entries in
$\boldsymbol{A}$ and two entries in $\boldsymbol{D}$. As a result, we need to
impose six identifiability constraints. We denote:

$\displaystyle\boldsymbol{A}=\begin{pmatrix}\zeta_{1}&\eta_{1}\\\
\zeta_{2}&\eta_{2}\end{pmatrix},\hskip
14.22636pt\boldsymbol{A}^{-1}=\dfrac{1}{\zeta_{1}\eta_{2}-\eta_{1}\zeta_{2}}\begin{pmatrix}\eta_{2}&-\eta_{1}\\\
-\zeta_{2}&\zeta_{1}\end{pmatrix},\hskip 7.11317pt\text{and}\hskip
7.11317pt\boldsymbol{D}=\begin{pmatrix}\delta_{1}\\\ \delta_{2}\end{pmatrix},$

with $\zeta_{i},\eta_{i},\delta_{i}\in\mathbb{R}$ for $i=1,2$ such that
$\zeta_{1}\eta_{2}-\eta_{1}\zeta_{2}\neq 0$. The transformations from Equation
(A.5) then result in:

$\displaystyle\tilde{B}^{(1)}_{x}$
$\displaystyle:=\zeta_{1}B_{x}^{(1)}+\zeta_{2}B_{x}^{(2)},\hskip
14.22636pt\tilde{K}^{(1)}_{t}:=\dfrac{1}{\zeta_{1}\eta_{2}-\eta_{1}\zeta_{2}}\left(\eta_{2}K_{t}^{(1)}-\eta_{1}K_{t}^{(2)}\right)+\delta_{1},$
(A.6) $\displaystyle\tilde{B}^{(2)}_{x}$
$\displaystyle:=\eta_{1}B_{x}^{(1)}+\eta_{2}B_{x}^{(2)},\hskip
14.22636pt\tilde{K}^{(2)}_{t}:=\dfrac{1}{\zeta_{1}\eta_{2}-\eta_{1}\zeta_{2}}\left(-\zeta_{2}K_{t}^{(1)}+\zeta_{1}K_{t}^{(2)}\right)+\delta_{2}.$

We then need to find the parameters
$\zeta_{i},\eta_{i},\delta_{i}\in\mathbb{R}$ such that the constraints in
Equations (A.2) and (A.4) are fulfilled.

##### Location constraints.

Equation (A.2) specifies the location constraint that applies to
$\tilde{K}_{t}^{(i)}$ for $i\in\\{1,2\\}$. Utilizing the transformation
outlined in Equation (A.6) in conjunction with the location constraint yields
the following result:

$\displaystyle\delta_{1}$
$\displaystyle=\dfrac{-1}{\zeta_{1}\eta_{2}-\eta_{1}\zeta_{2}}\left(\dfrac{\eta_{2}}{|\mathcal{T}|}\sum_{t\in\mathcal{T}}K_{t}^{(1)}-\dfrac{\eta_{1}}{|\mathcal{T}|}\sum_{t\in\mathcal{T}}K_{t}^{(2)}\right)$
(A.7) $\displaystyle\delta_{2}$
$\displaystyle=\dfrac{-1}{\zeta_{1}\eta_{2}-\eta_{1}\zeta_{2}}\left(-\dfrac{\zeta_{2}}{|\mathcal{T}|}\sum_{t\in\mathcal{T}}K_{t}^{(1)}+\dfrac{\zeta_{1}}{|\mathcal{T}|}\sum_{t\in\mathcal{T}}K_{t}^{(2)}\right).$

Next, we substitute the aforementioned expression into the transformed period
effects from Equation (A.6) and obtain:

$\displaystyle\tilde{K}_{t}^{(1)}$
$\displaystyle=\dfrac{1}{\zeta_{1}\eta_{2}-\eta_{1}\zeta_{2}}\left[\eta_{2}\accentset{\circ}{K}_{t}^{(1)}-\eta_{1}\accentset{\circ}{K}_{t}^{(2)}\right]$
(A.8) $\displaystyle\tilde{K}_{t}^{(2)}$
$\displaystyle=\dfrac{1}{\zeta_{1}\eta_{2}-\eta_{1}\zeta_{2}}\left[-\zeta_{2}\accentset{\circ}{K}_{t}^{(1)}+\zeta_{1}\accentset{\circ}{K}_{t}^{(2)}\right],$

where for $i\in\\{1,2\\}$:

$\accentset{\circ}{K}_{t}^{(i)}=K_{t}^{(i)}-\frac{1}{|\mathcal{T}|}\sum_{t\in\mathcal{T}}K_{t}^{(i)}.$

Note that, by construction, the transformed period effects in Equation (A.8)
satisfy the location constraints from Equation (A.2) for any $\eta_{1}$,
$\eta_{2}$, $\zeta_{1}$ and $\zeta_{2}$.

##### Orthogonal constraints.

We now put focus on the orthogonality constraints in Equation (A.4). Define
$\bar{\zeta}_{2}=\zeta_{2}/\zeta_{1}$ and introduce an identical notation for
$\bar{\eta}_{2}$. Using the notations in Equations (A.6) and (A.8), we obtain:

$\displaystyle\langle\tilde{\boldsymbol{B}}^{(1)},\tilde{\boldsymbol{B}}^{(2)}\rangle$
$\displaystyle=\zeta_{1}\eta_{1}\left\\{||\boldsymbol{B}^{(1)}||_{2}^{2}+\bar{\zeta}_{2}\bar{\eta}_{2}||\boldsymbol{B}^{(2)}||_{2}^{2}+(\bar{\zeta}_{2}+\bar{\eta}_{2})\langle\boldsymbol{B}^{(1)},\boldsymbol{B}^{(2)}\rangle\right\\}=0$
(A.9)
$\displaystyle\langle\tilde{\boldsymbol{K}}^{(1)},\tilde{\boldsymbol{K}}^{(2)}\rangle$
$\displaystyle=\dfrac{\zeta_{1}\eta_{1}}{(\zeta_{1}\eta_{2}-\eta_{1}\zeta_{2})^{2}}\left\\{-\bar{\zeta}_{2}\bar{\eta}_{2}||\accentset{\circ}{\boldsymbol{K}}^{(1)}||_{2}^{2}-||\accentset{\circ}{\boldsymbol{K}}^{(2)}||_{2}^{2}+(\bar{\zeta}_{2}+\bar{\eta}_{2})\langle\accentset{\circ}{\boldsymbol{K}}^{(1)},\accentset{\circ}{\boldsymbol{K}}^{(2)}\rangle\right\\}=0.$
(A.10)

We multiply Equation (A.9) with
$\langle\accentset{\circ}{\boldsymbol{K}}^{(1)},\accentset{\circ}{\boldsymbol{K}}^{(2)}\rangle$
and Equation (A.10) with
$\langle\boldsymbol{B}^{(1)},\boldsymbol{B}^{(2)}\rangle$ and then subtract
both equations to obtain:

$\displaystyle\overbrace{\big{(}\langle\accentset{\circ}{\boldsymbol{K}}^{(1)},\accentset{\circ}{\boldsymbol{K}}^{(2)}\rangle||\boldsymbol{B}^{(1)}||_{2}^{2}+\langle\boldsymbol{B}^{(1)},\boldsymbol{B}^{(2)}\rangle||\accentset{\circ}{\boldsymbol{K}}^{(2)}||_{2}^{2}\big{)}}^{:=\>\xi_{1}}\>+$
(A.11) $\displaystyle\hskip
142.26378pt\bar{\zeta}_{2}\bar{\eta}_{2}\underbrace{\big{(}\langle\accentset{\circ}{\boldsymbol{K}}^{(1)},\accentset{\circ}{\boldsymbol{K}}^{(2)}\rangle||\boldsymbol{B}^{(2)}||_{2}^{2}+\langle\boldsymbol{B}^{(1)},\boldsymbol{B}^{(2)}\rangle||\accentset{\circ}{\boldsymbol{K}}^{(1)}||_{2}^{2}\big{)}}_{:=\>\xi_{2}}=0.$

We proceed by substituting the above expression for $\bar{\zeta}_{2}$ into
Equation (A.10), resulting in a quadratic equation in $\bar{\eta}_{2}$:

$\displaystyle\overbrace{\left[\xi_{2}\langle\accentset{\circ}{\boldsymbol{K}}^{(1)},\accentset{\circ}{\boldsymbol{K}}^{(2)}\rangle\right]}^{:=\>a}\bar{\eta}_{2}^{2}+\overbrace{\left[\xi_{1}||\accentset{\circ}{\boldsymbol{K}}^{(1)}||_{2}^{2}-\xi_{2}||\accentset{\circ}{\boldsymbol{K}}^{(2)}||_{2}^{2}\right]}^{:=\>b}\bar{\eta}_{2}+\overbrace{\left[-\xi_{1}\langle\accentset{\circ}{\boldsymbol{K}}^{(1)},\accentset{\circ}{\boldsymbol{K}}^{(2)}\rangle\right]}^{:=\>c}=0.$
(A.12)

The corresponding discriminant equals:

$\displaystyle D$ $\displaystyle:=b^{2}-4ac$
$\displaystyle=\left[\xi_{1}||\accentset{\circ}{\boldsymbol{K}}^{(1)}||_{2}^{2}-\xi_{2}||\accentset{\circ}{\boldsymbol{K}}^{(2)}||_{2}^{2}\right]^{2}+4\xi_{1}\xi_{2}\langle\accentset{\circ}{\boldsymbol{K}}^{(1)},\accentset{\circ}{\boldsymbol{K}}^{(2)}\rangle^{2}.$

To prove that this discriminant is always positive, we consider two cases
based on the sign of $\xi_{1}\xi_{2}$. In the case $\xi_{1}\xi_{2}\geq 0$, it
is evident that $D\geq 0$. On the other hand, if $\xi_{1}\xi_{2}<0$, we employ
the Cauchy-Schwarz inequality and derive:

$\displaystyle D$
$\displaystyle\geq\left[\xi_{1}||\accentset{\circ}{\boldsymbol{K}}^{(1)}||_{2}^{2}-\xi_{2}||\accentset{\circ}{\boldsymbol{K}}^{(2)}||_{2}^{2}\right]^{2}+4\xi_{1}\xi_{2}||\accentset{\circ}{\boldsymbol{K}}^{(1)}||_{2}^{2}||\accentset{\circ}{\boldsymbol{K}}^{(2)}||_{2}^{2}$
$\displaystyle=\xi_{1}^{2}||\accentset{\circ}{\boldsymbol{K}}^{(1)}||_{2}^{4}+\xi_{2}^{2}||\accentset{\circ}{\boldsymbol{K}}^{(2)}||_{2}^{4}+2\xi_{1}\xi_{2}||\accentset{\circ}{\boldsymbol{K}}^{(1)}||_{2}^{2}||\accentset{\circ}{\boldsymbol{K}}^{(2)}||_{2}^{2}$
$\displaystyle=\left(\xi_{1}||\accentset{\circ}{\boldsymbol{K}}^{(1)}||_{2}^{2}+\xi_{2}||\accentset{\circ}{\boldsymbol{K}}^{(2)}||_{2}^{2}\right)^{2}$
$\displaystyle\geq 0.$

As a result, the quadratic equation in (A.12) has two real solutions:

$\displaystyle\bar{\eta}_{2}\in\left\\{\dfrac{-b-\sqrt{D}}{2a},\dfrac{-b+\sqrt{D}}{2a}\right\\}.$

By using the relation in Equation (A.11), we deduce the expression for
$\bar{\zeta}_{2}$:

$\displaystyle\bar{\zeta}_{2}=-\dfrac{\xi_{1}}{\xi_{2}\bar{\eta}_{2}}.$

Next, we prove that the solution sets of $\bar{\zeta}_{2}$ and
$\bar{\eta}_{2}$ coincide. For $\bar{\eta}_{2}=(-b-\sqrt{D})/2a$, we obtain:

$\displaystyle\bar{\zeta}_{2}$
$\displaystyle=\dfrac{2a\xi_{1}}{\xi_{2}\left(b+\sqrt{D}\right)}$
$\displaystyle=\dfrac{2\xi_{1}\xi_{2}\langle\accentset{\circ}{\boldsymbol{K}}^{(1)},\accentset{\circ}{\boldsymbol{K}}^{(2)}\rangle\left(b-\sqrt{D}\right)}{\xi_{2}\left(b^{2}-D^{2}\right)}$
$\displaystyle=\dfrac{-2\xi_{1}\xi_{2}\langle\accentset{\circ}{\boldsymbol{K}}^{(1)},\accentset{\circ}{\boldsymbol{K}}^{(2)}\rangle\left(b-\sqrt{D}\right)}{4\xi_{1}\xi_{2}^{2}\langle\accentset{\circ}{\boldsymbol{K}}^{(1)},\accentset{\circ}{\boldsymbol{K}}^{(2)}\rangle^{2}}$
$\displaystyle=\dfrac{-b+\sqrt{D}}{2\xi_{2}\langle\accentset{\circ}{\boldsymbol{K}}^{(1)},\accentset{\circ}{\boldsymbol{K}}^{(2)}\rangle}$
$\displaystyle=\dfrac{-b+\sqrt{D}}{2a}.$

Similarly, for $\bar{\eta}_{2}=(-b+\sqrt{D})/2a$, we can prove that
$\bar{\zeta}_{2}=(-b-\sqrt{D})/2a$ and vice versa.

##### Scaling constraints.

In the last step, we determine $\zeta_{1}$, $\eta_{1}$ such that the
constraints on the age effects in Equation (A.2) are fulfilled. We obtain:

$\displaystyle
1=\displaystyle\sum_{x\in\mathcal{X}}\left[\widetilde{B}_{x}^{(1)}\right]^{2}=\sum_{x\in\mathcal{X}}\left(\zeta_{1}B_{x}^{(1)}+\zeta_{1}\bar{\zeta}_{2}B_{x}^{(2)}\right)^{2}=\zeta_{1}^{2}\sum_{x\in\mathcal{X}}\left(B_{x}^{(1)}+\bar{\zeta}_{2}B_{x}^{(2)}\right)^{2}.$

As a result and using a similar reasoning for the constrained on the squared
$\widetilde{B}_{x}^{(2)}$, we deduce:

$\zeta_{1}=\dfrac{\text{sign}\left\\{\displaystyle\sum_{x\in\mathcal{X}}\left(B_{x}^{(1)}+\bar{\zeta}_{2}B_{x}^{(2)}\right)\right\\}}{\sqrt{\displaystyle\sum_{x\in\mathcal{X}}\left(B_{x}^{(1)}+\bar{\zeta}_{2}B_{x}^{(2)}\right)^{2}}},\hskip
14.22636pt\eta_{1}=\dfrac{\text{sign}\left\\{\displaystyle\sum_{x\in\mathcal{X}}\left(B_{x}^{(1)}+\bar{\eta}_{2}B_{x}^{(2)}\right)\right\\}}{\sqrt{\displaystyle\sum_{x\in\mathcal{X}}\left(B_{x}^{(1)}+\bar{\eta}_{2}B_{x}^{(2)}\right)^{2}}},$

where we include the $\text{sign}(\cdot)$ term to have a unique value for
$\zeta_{1}$ and $\zeta_{2}$. Note that the solution set of $\zeta_{1}$ and
$\eta_{1}$ again coincide. Since, e.g., $\zeta_{2}=\zeta_{1}\bar{\zeta}_{2}$,
we can easily obtain the (coinciding) solution sets of $\zeta_{2}$ and
$\eta_{2}$. These results also define the (coinciding) solution set of
$\delta_{1}$ and $\delta_{2}$ in Equation (A.7).

In conclusion, we have derived two different sets of parameter transformations
that lead to the same model fit and that satisfy the imposed location,
scaling, and orthogonality constraints. However, because of the coinciding
solution sets of $\eta_{1}$ and $\zeta_{1}$, $\eta_{2}$ and $\zeta_{2}$ and
$\delta_{1}$ and $\delta_{2}$, we obtain that the parameters
$\left\\{\tilde{B}_{x}^{(1)},\tilde{K}_{t}^{(1)},\tilde{B}_{x}^{(2)},\tilde{K}_{t}^{(2)}\right\\}$
in the first transformation coincides with the parameters
$\left\\{\tilde{B}_{x}^{(2)},\tilde{K}_{t}^{(2)},\tilde{B}_{x}^{(1)},\tilde{K}_{t}^{(1)}\right\\}$
in the second transformation. This correspondence does not pose a particular
issue as we can simply reassign the labels of the age and period effects. One
viable approach to ensure uniqueness entails selecting the parameter solution
such that the range of $\tilde{B}_{x}^{(1)}$ is larger than the one of
$\tilde{B}_{x}^{(2)}$.

##  Appendix B Indirect estimation of the mortality improvement model

We use the methodology proposed by Hunt and Villegas (2021) for indirect
estimation of the baseline mortality improvement model. Hereto, we first
transform the baseline mortality improvement model, as formulated in Equation
(3.1), into a baseline mortality model for the logarithmic force of mortality.
We iteratively calculate:

$\displaystyle\log\mu_{x,t}^{(c)}$
$\displaystyle=\log\mu_{x,t-1}^{(c)}+A_{x}+\sum_{i=1}^{m}B_{x}^{(i)}K_{t}^{(i)}+\sum_{j=1}^{l}\beta_{x}^{(j,c)}\kappa_{t}^{(j,c)}$
$\displaystyle=\log\mu_{x,t-2}^{(c)}+2A_{x}+\sum_{i=1}^{m}B_{x}^{(i)}\left(K_{t-1}^{(i)}+K_{t}^{(i)}\right)+\sum_{j=1}^{l}\beta_{x}^{(j,c)}\left(\kappa_{t-1}^{(j,c)}+\kappa_{t}^{(j,c)}\right)$
$\displaystyle=\ldots$
$\displaystyle=\log\mu_{x,t_{\min}}^{(c)}+\left(t-t_{\min}\right)A_{x}+\sum_{i=1}^{m}B_{x}^{(i)}\left(\sum_{\tau=t_{\min}+1}^{t}K_{\tau}^{(i)}\right)+\sum_{j=1}^{l}\beta_{x}^{(j,c)}\left(\sum_{\tau=t_{\min}+1}^{t}\kappa_{\tau}^{(j,c)}\right),$

for $x\in\mathcal{X}$ and $t\in\mathcal{T}$. We introduce the following
notation:

$\displaystyle
L_{t}^{(i)}=\displaystyle\sum_{\tau=t_{\min}+1}^{t}K_{\tau}^{(i)},\hskip
14.22636pt\lambda_{t}^{(j,c)}=\displaystyle\sum_{\tau=t_{\min}+1}^{t}\kappa_{\tau}^{(j,c)}.$
(B.1)

We then obtain an equivalent specification for the baseline mortality
improvement model in terms of the logarithmic force of mortality:

$\displaystyle\log\mu_{x,t}^{(c)}=\log\mu_{x,t_{\min}}^{(c)}+\left(t-t_{\min}\right)A_{x}+\displaystyle\sum_{i=1}^{m}B_{x}^{(i)}L_{t}^{(i)}+\displaystyle\sum_{j=1}^{l}\beta_{x}^{(j,c)}\lambda_{t}^{(j,c)},$
(B.2)

where $L_{t_{\min}}^{(i)}:=0$ for every $i\in\\{1,...,m\\}$ and
$\lambda_{t_{\min}}^{(j,c)}:=0$ for every $j\in\\{1,...,l\\}$.

Additionally, we need to impose identifiability constraints to make the
mortality model specification in Equation (B.2) fully identifiable. Remark
that the constraints formulated in Equations (3.2) and (3.3) serve as
identifiability constraints in the baseline mortality improvement model.
Consequently, we must translate these constraints in terms of the redefined
model specification outlined in Equation (B.2). Using Equation (B.1), we
translate the constraints in Equation (3.2) into:

$\displaystyle\sum_{x\>\in\>\mathcal{X}}\left(B_{x}^{(i)}\right)^{2}$
$\displaystyle=1,\hskip 14.22636ptL_{t_{\max}}^{(i)}=0,\hskip 28.45274pt$
$\displaystyle\text{for}\>i\in\\{1,\ldots,m\\}$ (B.3)
$\displaystyle\sum_{x\>\in\>\mathcal{X}}$
$\displaystyle\left(\beta_{x}^{(j,c)}\right)^{2}=1,$
$\displaystyle\text{for}\>j\in\\{1,\ldots,l\\}.$

The constraints in Equation (3.3) (for $m=2$) translate into:

$\displaystyle\displaystyle\sum_{x\>\in\>\mathcal{X}}B_{x}^{(1)}B_{x}^{(2)}=0,\hskip
14.22636pt\displaystyle\sum_{t=t_{\min}+1}^{t_{\max}}\left(L_{t}^{(1)}-L_{t-1}^{(1)}\right)\left(L_{t}^{(2)}-L_{t-1}^{(2)}\right)=0,$
(B.4)

with similar constraints on the country-specific parameters
$\beta_{x}^{(j,c)}$ and $\lambda_{t}^{(j,c)}$ for $j\in\\{1,2\\}$.

##  Appendix C Outlier detection using the Mahalanobis distance

We perform outlier detection using the Mahalanobis distance and apply a robust
Minimum Covariance Determinant (MCD) estimate of the location and scale. This
can be done in a three-step procedure:

  1. 1.

Compute the robust MCD estimate of location
$\hat{\boldsymbol{\mu}}^{\text{rob}}$ and scale
$\hat{\boldsymbol{\Sigma}}^{\text{rob}}$ for the $m$-dimensional time series
of remainder components
$\mathcal{R}=\left\\{\hat{\boldsymbol{R}}_{t}:=\left(\hat{R}_{t}^{(1)},...,\hat{R}_{t}^{(m)}\right)\mid
t\in\mathcal{T}\right\\}$, i.e. we use the CovMcd function from the rrcov
package in R.

  2. 2.

Compute the robust distances for each observation $\hat{\boldsymbol{R}}_{t}$:

$\displaystyle
d_{t}=\sqrt{\left(\hat{\boldsymbol{R}}_{t}-\hat{\boldsymbol{\mu}}^{\text{rob}}\right)^{T}\left(\hat{\boldsymbol{\Sigma}}^{\text{rob}}\right)^{-1}\left(\hat{\boldsymbol{R}}_{t}-\hat{\boldsymbol{\mu}}^{\text{rob}}\right)}.$

  3. 3.

We identify year $t$ as an outlier whenever the robust distance $d_{t}$
exceeds the square root of the $99\%$ quantile of the chi-square distribution
with $m$ degrees of freedom, i.e.
$\sqrt{\scalebox{1.5}{$\chi$}_{m,0.99}^{2}}$.

##  Appendix D A weighted maximum likelihood estimation method for
projecting common and country-specific period effects

We provide more details on the approach used for estimating the multivariate
time series model for the joint vector of common and country-specific period
effects, as formulated in Equation (6.7). We implement a weighted variant of
maximum likelihood estimation. Hereto, we minimize the following objective
function:

$\displaystyle-\displaystyle\sum_{t=1851}^{t_{\max}}\omega_{t}\log\left(\phi\left(\boldsymbol{\mathcal{K}}_{t}\mid\boldsymbol{c},\boldsymbol{\Sigma}_{w}\right)\right),$
(D.1)

where $\boldsymbol{\mathcal{K}}_{t}$ is the four-dimensional vector of period
effects, $w_{t}=\gamma^{t_{\max}-t}$ are geometrically decaying weights and
$\gamma$ is the rate of decay.292929It should be noted that when the value of
the parameter $\gamma$ is set equal to one, the resulting procedure is
equivalent to traditional maximum likelihood estimation. Furthermore,
$\phi(\>\cdot\mid\boldsymbol{c},\boldsymbol{\Sigma}_{w})$ denotes the four-
dimensional Gaussian density function with mean $\boldsymbol{c}$ and
covariance matrix $\boldsymbol{\Sigma}_{w}$. The decay parameter $\gamma$
needs to be selected properly since it will have a major impact on the final
projections. We follow an approach similar to Mittnik and Paolella (2000) and
select the parameter $\gamma$ that maximizes the average of the log-likelihood
values associated with the one-step ahead projections since the year 1900
onwards,303030We start the one-step ahead projections at the year 1901, to
have sufficient data available to calibrate the time series models, i.e. at
least 50 years. i.e. we maximize:

$\displaystyle\hat{\gamma}=\operatorname*{arg\,max}\limits_{\gamma}\frac{1}{121}\sum_{t=1900}^{2020}\log\left(\hat{\phi}_{t+1\mid
t}\left(\boldsymbol{\mathcal{K}}_{t+1}\mid\gamma,\hat{\boldsymbol{c}},\hat{\boldsymbol{\Sigma}}_{w},\boldsymbol{\mathcal{K}}_{t},\ldots,\boldsymbol{\mathcal{K}}_{1851}\right)\right).$
(D.2)

Here, we calculate the one-step ahead projection by estimating the time series
models up to the year $t$, followed by projecting the calibrated model to year
$t+1$. Hence, the time series parameters $\hat{\boldsymbol{c}}$ and
$\hat{\boldsymbol{\Sigma}}_{w}$ in Equation (D.2) are re-estimated for each
value of $t$. The optimally selected value for $\gamma$ in the case study of
Section 6 equals $0.943$.

##  Appendix E Calibration and projection the regime-switching model

Our proposed calibration strategy builds upon prior work of Hainaut (2012).
However, we modify his methodology to account for the age-specific mortality
shock effect. Furthermore, to ensure that the Markov chain remains in the high
volatility regime for at least two subsequent years, we introduce a memory
state. This is necessary because a mortality shock in a given year $t$ leads
to outlying values in the mortality improvement rates for both year $t$ and
year $t+1$. In Section E.2, we discuss projections from our regime-switching
model. Specifically, we generate artificial mortality shocks in a manner that
aims to reflect real-world scenarios.

###  E.1 Calibration

####  E.1.1 Setting up the log-likelihood

The log-likelihood of the regime-switching model equals:

$\displaystyle l(\Theta)$ $\displaystyle=\log
f\left(\boldsymbol{z}_{t_{\min}},\boldsymbol{z}_{t_{\min+1}},...,\boldsymbol{z}_{t_{\max}}\mid\Theta\right)$
(E.1) $\displaystyle=\sum_{t\in\mathcal{T}}\log
f\left(\boldsymbol{z}_{t}\mid\boldsymbol{z}_{t-1},...,\boldsymbol{z}_{t_{\min}},\Theta\right),$

with $\mathcal{T}:=\left\\{t_{\min},t_{\min}+1,\ldots,t_{\max}\right\\}$ the