1. Introduction
In many applications, there is a need for reliable mortality models. For example, they are an essential ingredient of population projections, which yield important information for planning the future of society. For pension funds and life insurers, accurate mortality modelling is part of their business model.
The last century has seen a tremendous decrease of mortality rates in developed countries, calling attention to the risk of underestimating average future longevity, which has actually happened in the official estimates of many of these countries (Michaelson & Mulholland, Reference Michaelson and Mulholland2014). This risk has been termed longevity risk and particularly concerns all types of annuity providers as their liabilities depend on mortality development. Many efforts such as improvements in modelling and the design of innovative capital market hedging instruments have been undertaken to cope with longevity risk.
Recently, a threat to the continued longevity improvements has emerged in the shape of the COVID-19 pandemic, which has cost over 4.6 million lives worldwide as of September 2021 (Dong et al., Reference Dong, Du and Gardner2021). SARS-CoV-2 is not the first virus to cause a pandemic, and it certainly will not be the last. A study by Woolhouse & Gaunt (Reference Woolhouse and Gaunt2007) finds that over 6% of all known human pathogens were first reported in humans after 1980, and many of them were globally spreading viruses.
Any mortality modeller must address the question whether allowing for pandemics, or more generally for mortality shocks induced not only by pandemics but also other events such as natural catastrophes or major terror attacks, is necessary for their specific application. In contrast to annuity providers, writers of life insurance policies are exposed to mortality risk, i.e., the risk of underestimating future mortality rates. This risk might be enlarged by not taking into account the possibility of mortality jumps. Therefore, as Cox et al. (Reference Cox, Lin and Pedersen2010) write, “including pandemic effects is an important issue in modeling mortality for life insurance liabilities”.
Accounting for severe mortality shocks is in fact required from life insurers by Solvency II regulations, according to which the available mortality risk capital must be sufficient to cover for losses resulting from an instantaneous permanent increase of mortality rates by 15% (European Union, 2020). While it seems unlikely that COVID-19 will cause such an extreme long-term change of mortality rates, it is plausible that even transitory mortality jumps of this size can have non-negligible financial consequences.
An empirical study on US life insurance data by Individual Life COVID-19 Project Work Group (2021) finds that life insurance claim counts have increased by 4.6% for females and 5.5% for males in the first half of 2020. The appearance of excess claims differs by age and heavily by geographical region. In particular, the New York/New Jersey region “had a terrible April”: all-cause life claim counts were 230% of the previous five years’ average. Similar analyses have been performed in other countries. For example, Continuous Mortality Investigation (2020) present empirical evidence for weekly excess mortality of up to 144% in England and Wales along with some recommendations on whether and how to include 2020 data in the mortality modelling and forecasting process. They particularly mention the necessity to adjust stochastic mortality models but do not provide a deeper analysis.
In the following, we investigate whether COVID-19 had a notable impact on mortality rates in different regions by analysing mortality data up to and including 2020 for several European countries. We find that the influence of COVID-19 on mortality in 2020 is clearly recognisable both on a weekly and on a yearly scale according to different measures such as excess death ratios and age-standardised death rates. COVID-19 has caused a (temporary) setback of the mortality development in the considered countries by up to 12 years. In terms of improvement rates, 2020 is among the 10 worst years for all countries in our considered data set, which dates back to as far as 1900 in many cases. Therefore, it seems advisable to take COVID-19 into account in mortality models and for valuations, reserving decisions or solvency capital calculations derived from such models by (re)insurance companies.
In order to quantify the impact of the 2020 mortality shock, we calibrate the widely used Lee–Carter (LC) model by Lee & Carter (Reference Lee and Carter1992) on real data including COVID-19 deaths and compare its parameters, forecasts and the implied present values of annuities and life insurance policies to a counterfactual scenario in which 2020 mortality continues to develop according to the 2019 trend. As expected, we find that forecasts based on the real data drift away from forecasts based on the hypothetical situation without COVID-19 with increasing forecasting horizon. Consequently, annuity values decrease by up to 9% and term assurance values increase by up to 29% in the numerical example we consider. This is due to the fact that we use 2020 as the forecast jump-off year, whose mortality strongly influences the LC forecasts under the usual random walk with drift forecasting method. However, with another numerical example, we show that a mortality jump at any time point in the calibration data set of an LC model has a significant impact on the uncertainty related to its forecasts and should therefore not be ignored but explicitly dealt with by the modeller. To check whether our findings pertain to other mortality models, we consider the Cairns–Blake–Dowd (CBD) model by Cairns et al. (Reference Cairns, Blake and Dowd2006) as a robustness check. We obtain qualitatively similar results.
Additionally, we provide a theoretical analysis regarding the influence of a mortality shock during the calibration period on the parameters and forecasts of the LC model. To facilitate understanding of our results, we make an assumption on the shape of this shock. In the demographic literature, evidence has been found that the age dependence of COVID-19-related mortality is similar across a broad range of countries and often resembles the age dependence of all-cause mortality (Goldstein & Lee, Reference Goldstein and Lee2020). There is some variation in COVID-19 mortality age patterns between countries, which is likely due to country-specific factors such as health care standards, but Sasson (Reference Sasson2021) empirically confirms that COVID-19 mortality in many countries increases approximately log-linear with age between ages 25 and 84. Therefore, we will base our theoretical analysis on the assumption that COVID-19 (approximately) induces a log-parallel shift of death rates, i.e., it makes mortality rates increase by an age-independent multiplicative factor (see also Milevsky, Reference Milevsky2020).
To summarise, we illustrate the influence of COVID-19 on mortality rates and quantify its potential impact on insurance valuation applications. We find this impact to be quite high for point forecasts when using 2020 as the jump-off year and we illustrate that the impact on interval forecasts will continue to be high, even if mortality levels were to normalise quickly. Speaking in statistical terms, a mortality jump might lead to a substantial violation of the usual normal distribution assumption made for the period effect increments of the LC model. This implies that the random walk with drift, which is the most commonly used time series model for obtaining forecasts in the LC framework, might lead to invalid conclusions and unreliable forecasts when based on data which include a mortality shock. This observation leads the way to several approaches from mathematical finance to deal with the situation, of which we provide an overview as a starting point for further research.
We proceed as follows: In Section 2, the different sources of data used in our study along with some necessary preprocessing steps are described. This is followed by an explanation of the methods used for analysing mortality data and the LC model as well as theoretical considerations regarding the impact of a mortality jump during the model calibration period in Section 3. The results of our empirical analysis are given in Section 4. In Section 5, we present several approaches from the literature regarding the adjustment of mortality data or models to account for (the possibility of) mortality shocks. Section 6 concludes. Details on the LC model, the proofs of all theoretical findings and some additional empirical results are presented in the appendix.
2. Data
We analyse mortality data which depend on age
$x\in\mathcal{X} = \{x_1, \dots, x_A\}$
, where A denotes the number of ages or age groups. For the numerical results in Section 4, the available data are given in 5-year age groups starting from age 35 so that
$x_1 = 35$
,
$x_2 = 40, \dots, x_A = 90$
and
$A = 12$
. Note that we denote the age groups by their lower bound for notational simplicity. Explicitly, we consider the age intervals
$35-39, 40-44, \dots, 85-89$
and the open-ended age group
$90+$
. The data also depend on the calendar year, which we denote by
$t \in \mathcal{T} = \{t_1, \dots, t_Y\}$
, where Y is the number of available years and differs by population as detailed below.
We consider the mortality experience of multiple populations which we index by
$i\in\mathcal{P} \;:\!=\; \{1,\dots, P\}$
, where P denotes the number of populations. In our application, a population is a combination of country and sex (female, male, total). We restrict our analysis to Austria (1947–2020), Belgium, France, Germany (1956–2020, aggregating West and East Germany before 1990), Italy, Poland (1958–2020), Spain (1908–2020), Sweden and Switzerland. Where not stated differently, data are available for the period 1900–2020. This means that we consider
$P = 27$
populations in total.
Our main quantity of interest is the death rate, which is the proportion
\begin{equation} m_{x,\;t}^i \;:\!=\; \frac{D_{x,\;t}^i}{E_{x,\;t}^i}\end{equation}
of the death count
$D_{x,\;t}^i$
for population i, year t and age x to the corresponding total number of person-years at risk, the exposure
$E_{x,\;t}^i$
. We also consider data on a weekly scale and write
$m_{x,\;t,\;w}^i$
,
$D_{x,\;t,\;w}^i$
and
$E_{x,\;t,\;w}^i$
with
$w\in\{1,\dots,52\}$
for weekly death rates, death counts and exposures.
We get observations of these quantities from the Human Mortality Database (2021, HMD) and Short Term Mortality Fluctuations (2021, STMF). They provide historical yearly and recent weekly mortality data, respectively, of several countries in a unified format. Note that there can be substantial delay in the reporting of fully accurate vital statistics, which means the most recent STMF data should be considered preliminary. However, as a consequence of the increased public interest in mortality data during the COVID-19 pandemic, many statistical offices have made significant efforts to provide high-quality data in a more timely manner. We expect the effects of future data updates to be negligible for our conclusions
Before we use the data in our analysis, we apply some preprocessing steps. We aggregate death counts and exposures where necessary to obtain the 5-year age groups as described above. In the STMF data for Sweden, there are some deaths for which age or week of occurrence are unknown. We redistribute deaths at unknown ages according to the observed distribution of deaths at known ages, and we also redistribute deaths in unknown weeks according to the observed distribution of deaths in known weeks within each year. For details on this standard approach, we refer to the HMD methods protocol (Wilmoth et al., Reference Wilmoth, Andreev, Jdanov, Glei and Riffe2021, section 4.1).
Furthermore, STMF data are available on a weekly scale, but for some purposes we need recent yearly mortality data. For such purposes, we group the weekly death counts in the STMF data by age, population and year and aggregate them to obtain yearly death counts. These yearly death counts are then appended to the HMD data where the most recent years (2018–2020, 2019–2020 or only 2020) are not available yet. We expect this to give a fairly accurate impression on how recent trends compare to historical observations.
Finally, to calculate death rates from death counts, we also need the corresponding exposure, which is not available in the STMF input data. As exposures are notoriously hard to estimate, we refrain from sophisticated methods and perform a simple age- and population-specific linear extrapolation of the last five known exposures from the HMD (similarly to Leavitt, Reference Leavitt2021). This will not be completely accurate but should be sufficient to get a good qualitative understanding of the recent mortality developments. Setting
$E_{x,\;t,\;w}^i \;:\!=\; \frac{E_{x,\;t}^i}{52}$
for all the weeks
$w=1,\dots,52$
in a year
$t\leq 2020$
allows us to calculate weekly death rates as well.
3. Methodology and Theoretical Results
3.1. Mortality statistics
To quantify the impact of COVID-19 on mortality, we calculate several measures derived from the raw death counts or death rates. A term which has often been used in the context of the pandemic is excess mortality, relating to the question by how much mortality during the COVID-19 pandemic deviates from previous expectations. We prefer to investigate all-cause excess mortality as opposed to directly analysing data on COVID-19-related deaths because these might give an incomplete account for multiple reasons (such as cause of death misclassifications and collateral effects of the pandemic on other causes of death, see Aburto et al., Reference Aburto, Kashyap, Schöley, Angus, Ermisch, Mills and Dowd2021; Institute for Health Metrics and Evaluation 2021b for further details). There are several ways to measure excess mortality (Aron & Muellbauer, Reference Aron and Muellbauer2020). We will look at weekly excess death ratios, defined by
\begin{equation} p_{x,\;t,\;w}^i \;:\!=\; \frac{D_{x,\;t,\;w}^i - \hat{D}_{x,\;t,\;w}^i}{\hat{D}_{x,\;t,\;w}^i},\end{equation}
where
$\hat{D}_{x,\;t,\;w}^i$
denotes our expectation under “normal” circumstances for the death counts in week w of year
$t\in\mathcal{T}$
. An often used approach is to define this as the average of the
$N_P$
previous years (Eurostat, 2021; Destatis, 2021), i.e.,
\begin{equation} \hat{D}_{x,\;t,\;w}^i \;:\!=\; \frac{1}{N_P} \sum\limits_{l=1}^{N_P} D_{x,\;t-l, w}^i,\end{equation}
where we use
$N_P \;:\!=\; 4$
. This average does not necessarily coincide with our expectations, as one might rather expect mortality to decrease over time and as it does not take into account changes in population size and structure, but we deem (3) sufficient to get a qualitative impression.
We also look at mortality rates averaged over ages. However, in order to make meaningful comparisons over different populations with potentially quite different age structures, we apply a standardisation procedure yielding age-standardised death rates
\begin{equation} m_{S,\;t}^{i} \;:\!=\; \sum\limits_{x=x_1}^{x_A} \omega_x m_{x,\;t}^i,\end{equation}
where the weights
$\omega_x\in[0,1]$
are defined by
\begin{equation} \omega_x \;:\!=\; \frac{E_x^S}{\sum\limits_{x=x_1}^{x_A} E_x^S}.\end{equation}
Here,
$E_x^S$
denotes exposures of a synthetic population, in this case the European Standard Population (European Commission & Eurostat, 2013). This makes age-standardised death rates independent of the underlying age structure and therefore facilitates comparability over populations (Keyfitz & Caswell, Reference Keyfitz and Caswell2005, Chapter 4.1).
To achieve comparability over time as well, a standard approach is to consider changes in mortality rates instead of the mortality rates themselves. We implement this by calculating annual improvement rates
\begin{equation} I_{x,\;t}^i \;:\!=\; \frac{m_{x,\;t-1}^i - m_{x,\;t}^i}{m_{x,\;t-1}^i},\end{equation}
which are positive if mortality decreases from
$t-1$
to t and negative if it increases. We again perform age standardisation via
\begin{equation} I_{S,\;t}^i \;:\!=\; \sum\limits_{x=x_1}^{x_A} \omega_x I_{x,\;t}^i\end{equation}
with the weights as in (5)Footnote 1 .
3.2. Mortality forecasts and insurance contracts
The mortality model by Lee & Carter (Reference Lee and Carter1992) has become a standard both in academic work and applications. Therefore, we mainly focus our analysis on this model in order to evaluate the influence of COVID-19 on mortality modelling. For the reader’s convenience, we provide some details on model calibration and on obtaining point and interval forecasts in Appendix A.
Of course, the LC model has some limitations, for example, it does not take into account cohort effects, which are certainly relevant for some of the populations we consider. The following evaluation in Section 4 could easily be repeated for the LC model with cohort effects (Renshaw & Haberman, Reference Renshaw and Haberman2006) or, in fact, any other stochastic mortality model. We do not expect that this would change our main qualitative conclusions, which will be confirmed by a robustness check with the Cairns–Blake–Dowd model in Section 4.5.
In order to quantify the impact of COVID-19 on mortality forecasts, we calibrate LC models on different data sets. An overview is given in Table 1.
Overview of the different LC models we use and their calibration periods.
In particular, we will present two comparative analyses, the first of which is between
-
• real_20 (data from 1991 to 2020, where 2020 data are observed mortality rates), and
-
• LC_20 (data from 1991 to 2020, where 2020 data are the central forecasts of an auxiliary LC model trained on the years 1991 to 2019).
In this case, using the notation from Section 2, we have
$\mathcal{T} = \{1991,\dots,2020\}$
and
$Y=30$
for all populations. Analysing these two models, we get an impression of the behaviour of a model including COVID-19 mortality compared to a model trained in a world where 2020 mortality behaves exactly as the LC model would have expected in 2019. This illustrates the consequences of a mortality shock in the jump-off year of the forecast. In other words, COVID-19 has influenced mortality in 2020 and we calibrate our mortality model on data up to and including 2020 and make forecasts from this year on. It is clear from the forecasting method, see (A.4), that mortality in the last year
$t_Y$
has a decisive impact on point forecasts.
If mortality returned to normal levels in a future year
$t > 2020$
, central forecasts based on data including t would not be influenced much by the mortality jump in 2020. However, whether the jump is present or not would certainly make a difference for the prediction intervals (A.5) because it increases the estimated volatility of the period effect time series. In order to study these two effects with a numerical example, we will make a second comparison of two LC models calibrated on slightly different data sets:
-
• real_20/LC_21 (data from 1992 to 2021, where 2020 data are observed mortality rates and 2021 data are the central forecasts of an auxiliary LC model trained on the years 1991 to 2019), and
-
• LC_20/21 (data from 1992 to 2021, where both 2020 and 2021 rates are obtained from the auxiliary LC model).
In this case, using the notation from Section 2, we have
$\mathcal{T} = \{1992,\dots,2021\}$
and
$Y=30$
for all populations. The aim is to compare the real development including COVID-19 and an assumed subsequent return to the expected levels as of 2019 on the one hand with a world where COVID-19 did not occur and death rates consistently behave according to the LC model on the other hand. Note that we do not claim that 2021 mortality is going to return to LC levels. We just consider it as a scenario to illustrate how the LC forecasts are influenced by a mortality shock which happens within the calibration period but not in the jump-off year. We expect the obtained results to be qualitatively similar to the situation where a return to normal levels only occurs some time after 2021.
We measure how COVID-19 mortality could influence insurance contract values by calculating present values of annuities and term life insurance policies under the different LC models. We assume a constant yearly discount factor v. This is a simplification allowing us to focus solely on mortality shocks. First, we consider a temporary annuity-immediate issued to a life aged x at the start of year t which runs for n years and pays an amount of 1 at the end of each year in which the annuitant is alive. Denoting the s-year cohort survival probability by
$_{s}p_{x,\;t}$
with the special case
$p_{x,\;t} \;:\!=\; {}_{1}p_{x,\;t}$
, the present value of such an annuity is given by
\begin{equation} a_{x:\overline{n}\kern-0.5pt\raise.5pt\hbox{${\scriptstyle{\mid}}$}}(t) \;:\!=\; \sum\limits_{s=1}^n v^s {}_{s}p_{x,\;t} =\sum\limits_{s=1}^n v^s \prod\limits_{j=0}^{s-1} p_{x+j,\;t+j} \approx \sum\limits_{s=1}^n v^s \exp\left(-\sum\limits_{j=0}^{s-1} m_{x+j, t+j}\right).\end{equation}
Furthermore, we consider a term assurance issued to a life aged x at the start of year t which runs for n years and pays an amount of 1 at the end of the year if the life has died within this year. Its present value is given by
\begin{equation} \begin{aligned} A_{x:\overline{n}\kern-0.5pt\raise.5pt\hbox{${\scriptstyle{\mid}}$}}(t) &\;:\!=\; \sum\limits_{s=0}^{n-1} v^{s+1} {}_{s}p_{x,\;t} \left(1 - p_{x+s,\;t+s}\right) \\[5pt] &\approx \sum\limits_{s=0}^{n-1} v^{s+1} \exp\left(-\sum\limits_{j=0}^{s-1} m_{x+j, t+j}\right) \left(1 - \exp\left( -m_{x+s,\;t+s} \right)\right). \end{aligned}\end{equation}
In the above equations, we have used the standard approximation
$p_{x,\;t} \approx \exp\left(-m_{x,\;t}\right)$
(Pitacco et al., Reference Pitacco, Denuit, Haberman and Olivieri2008, Chapter 2.3).
We calculate intervals for annuity and life insurance values simply by inserting the prediction interval bounds for the death rates from (A.5) into the valuation formulae (8) and (9). For the annuity, the upper bound is obtained by inserting the lower bounds of the death rates, and vice versa, due to the monotonically decreasing dependence of annuity values on death rates.
3.3. The impact of a log-parallel mortality shift on model parameters and forecasts
We close this section with a theoretical analysis of the question how a mortality shock influences the parameters and forecasts of an LC model. Generally, being able to estimate by how much a mortality model changes in response to the arrival of new calibration data is highly relevant in applications. For instance, under the Solvency II framework, solvency capital requirements are calculated as the 1-year value-at-risk of the basic own funds of an insurance company. When estimating the distribution of own funds in 1 year’s time, the fact that new mortality data will become available during this year has to be taken into account for valuing mortality-related liabilities. New data would lead to a recalibration of the applied mortality model with the risk of a subsequent significant change in its forecasts and, therefore, the valuation of liabilities. This risk has been termed recalibration risk and demonstrated to be material for the LC model (Cairns, Reference Cairns2013).
Here, we exemplarily focus on the special case that mortality rates at all ages increase by a common multiplicative factor, which is equivalent to a parallel shift of logarithmic mortality rates (“log-parallel shift”). As detailed in Section 1, this has been observed to approximately be the case for the COVID-19 mortality shock in several populations between ages 25 and 89. Furthermore, it is also the scenario used in the Solvency II standard formula for mortality risk quantification.
In the following, we suppress the dependence of the model parameters on the population i for better readability. We write
$m_{x,\;t}^B$
for the baseline mortality within a population. We assume
$m_{x,\;t} = m_{x,\;t}^B$
for all
$x\in\mathcal{X}$
and
$t < t_Y$
, whereas baseline mortality is not observed during the extreme event which we assume takes place in the final year
$t_Y$
and increases mortality by
$c \cdot m_{x,\;t_Y}^B$
for some
$c > 0$
, so that
$m_{x,\;t_Y} = (1+c) m_{x,\;t_Y}^B$
. Note that this implies
\begin{equation} \frac{m_{x,\;t_Y} - m_{x,\;t_Y}^B}{m_{x,\;t_Y}^B} = c,\end{equation}
so that c is the age-independent excess mortality ratio with respect to baseline mortality. As explained above, we assume c to be independent of age, which makes the following analysis easier to grasp but is not strictly necessary, see Remark 4 below.
We consider an LC model for the observed rates
\begin{equation} \log m_{x,\;t} = \alpha_x + \beta_x \kappa_t + \varepsilon_{x,\;t},\end{equation}
and an LC model for the baseline rates
\begin{equation} \log m_{x,\;t}^B = \alpha_x^B + \beta_x^B \kappa_t^B + \varepsilon_{x,\;t}^B.\end{equation}
We write
$\Delta \varepsilon_{x,\;t} \;:\!=\; \varepsilon_{x,\;t} - \varepsilon_{x,\;t}^B$
and assume that c is small enough to imply
\begin{equation} \Delta \varepsilon_{x,\;t} \approx 0 \text{ for all } x\in\mathcal{X}, t\in\mathcal{T}.\end{equation}
A more explicit formulation of (13) is that the exceedance probability
$\mathbb{P}\left(\left| \Delta\varepsilon_{x,\;t} \right| > \delta \right)$
is required to be low for some small
$\delta > 0$
, which is achieved if the variance
$\sigma_{\Delta\varepsilon}^2$
is sufficiently close to 0 by Chebyshev’s inequality. In this sense, assumption (13) and the following results could be made more precise, but we refrain from doing so for clarity of exposition. Proofs of the following statements are given in Appendix C.
Proposition 1 (Parameter changes induced by a log-parallel shift). Under the identifiability constraints
\begin{equation} \sum\limits_{x=x_1}^{x_A} \beta_x = 1 \text{ and } \kappa_{t_1} = 0 \end{equation}
and assumption (13), we get the following relationships between the parameters of the LC models (11) and (12):
\begin{equation} \alpha_x \approx \alpha_x^B, \end{equation}
\begin{equation} \beta_x \approx \beta_x^B, \end{equation}
\begin{equation} \kappa_t \approx \kappa_t^B + A\log\left(1+c\right) {\unicode[Times]{x1D7D9}}_{t=t_Y}, \end{equation}
where
$\approx$
denotes equality between both sides up to (terms of)
$\Delta\varepsilon_{x,\;t}$
.
Proposition 1 confirms the intuition that a mortality shock in a particular year t should mainly lead to a significant change (an increase) of the period effect
$\kappa_t$
belonging to that year. The details of the result depend on the choice of identifiability constraints, but it is easily transferred to other constraints and its qualitative statement holds for all constraints which are usually applied in the literature.
As the parameters of a model determine its forecasts, we can analyse these as well.
Corollary 1. (Forecast changes induced by a log-parallel shift). Under the conditions of Proposition 1 and using the random walk with drift forecasting approach, we get the following relationship between the h-step central forecasts of the LC models (11) and (12):
\begin{equation} \hat{m}_{x,\;t_Y+h} \approx \hat{m}_{x,\;t_Y+h}^B \cdot \exp\left(\hat{\beta}_{x}^B \cdot A \log\left(1+c\right) \left(1 + \frac{h}{Y-1}\right) \right). \end{equation}
Remark 1 (Interpretation of Corollary 1). Corollary 1 shows the influence of different parameters on the LC forecasts when a shock occurs in the jump-off year. A higher age effect
$\hat{\beta}_{x}^B$
leads to an increased sensitivity of the model death rate to temporal mortality development and, thus, to a larger increase in response to a shock. Unsurprisingly, the ratio of
$\hat{m}_{x,\;t_Y+h}$
and
$\hat{m}_{x,\;t_Y+h}^B$
also increases with the excess death ratio c, the number of available ages A and the length of the forecasting horizon h. On the other hand, the influence of a shock on the forecast can be diminished by increasing the number Y of training years. However, even with a very long history of training data mortality forecasts will increase approximately by a factor of
$\exp\left(\hat{\beta}_x^B \cdot A \log\left(1+c\right)\right)$
as a result of the shock in the jump-off year.
To give a numerical example, the 1-step forecast at age group
$x = 60$
with age effect
$\hat{\beta}_x^B = 0.075$
increases by
$9.3\%$
(compared to the situation without shock) if there are
$A = 12$
age groups,
$Y = 30$
years of training data and a
$c = 10\%$
shock occurs in the jump-off year
$t_Y$
. Ceteris paribus, the increase in the 20-step forecast is
$15.6\%$
.
Remark 2 (Generalisation to prediction intervals and present values). Corollary 1 can be generalised to describe the shock-induced change of the LC random walk variance estimator (A.6). This would allow us to make statements about the change in uncertainty related to the point forecasts (A.3). Similarly, building on Corollary 1, the impact of a mortality shock on annuity and life insurance values as defined in (8) and (9) can be quantified. We refrain from deriving explicit formulae because these would be rather complex and give little additional insight.
Remark 3 (Generalisation to other shock times). So far, we have assumed that a mortality shock occurs at time
$t_Y$
, i.e., the jump-off year of the mortality forecasts. Of course, the situation where a mortality shock occurs sometime between
$t_1$
and
$t_Y$
is of interest as well. Proposition 1, Corollary 1 and Remark 2 are easily generalised to a mortality jump occurring at any time
$t\in\{t_2,\dots,t_{Y-1}\}$
.
In particular, if the shock occurs at
$t<t_Y$
, random walk forecasts as in Corollary 1 should change very little compared to the situation without a jump because they depend mainly on the drift
$\hat{\mu}$
, which is determined by the marginal period effects
$\hat{\kappa}_{t_1}$
and
$\hat{\kappa}_{t_Y}$
. However, the uncertainty in the forecasts, which additionally depends on the variance of the random walk, would increase. Both these effects were observed by Lee & Carter (Reference Lee and Carter1992) when calibrating their model on data including the 1918 pandemic.
Remark 4 (Generalisation to non-parallel shifts). Despite the empirical evidence that COVID-19 might have caused an approximate log-parallel shift of the death rates, this is obviously not the case for every extreme mortality event. The above results can be generalised in this direction as well by allowing c in (10) to depend on age x, replacing c with
$c_x$
and repeating the calculations performed in the proofs.
Weekly, country-specific age-standardised death rates
$m_{S,\;t,\;w}^i$
as defined in (4) for the years 2019 and 2020.
Weekly, country-specific excess death ratios
$p_{x,2020,w}^i$
as defined in (2) for the year 2020 and age groups 40–44, 50–54, 60–64, 70–74, 80–84 and 90+. Values above the zero line (blue, dash) indicate excess mortality.
As a further illustration, we provide the analogue of Proposition 1 for a variant of the CBD model,
\begin{equation} \log m_{x,\;t} = \kappa_t^{(1)} + (x-\bar{x}) \kappa_t^{(2)} + \varepsilon_{x,\;t},\end{equation}
where
$\bar{x} \;:\!=\; \frac{1}{A} \sum\limits_{x=x_1}^{x_A} x$
and
$\varepsilon_{x,\;t} \sim \mathcal{N}\left(0, \sigma^2\right)$
i.i.d.
Proposition 2 (Parameter changes of the CBD model induced by a log-parallel shift). For a CBD model (19) of observed mortality
$m_{x,\;t}$
and another CBD model of baseline mortality
$m_{x,\;t}^B$
and given an assumption similar to (13), it holds that
\begin{align} \kappa_t^{(1)} \approx \kappa_t^{(1),B} + \log\left(1+c\right){\unicode[Times]{x1D7D9}}_{t=t_Y}, \end{align}
\begin{align} \kappa_t^{(2)} \approx \kappa_t^{(2),B}. \end{align}
4. Empirical Results
4.1. Mortality development in 2020
Figure 1 shows weekly age-standardised death rates
$m_{S,\;t,\;w}^i$
during the years 2019 and 2020. Dramatic increases are visible in some countries at the time when COVID-19 was emerging in Europe and in the winter months at the end of 2020.
In order to quantify exactly how dramatic these increases were, we compare the development of 2020 to the average of the four preceding years by calculating excess death ratios
$p_{x,2020,w}^i$
, which are displayed in Figure 2. For a better overview, we only show some of the age groups. Furthermore, we focus on the total populations, although men have a significantly higher risk of dying from a COVID-19 infection compared to women (Jin et al., Reference Jin, Bai, He, Wu, Liu, Han, Liu and Yang2020). Excess mortality is clearly visible at higher ages for all populations, with Belgium, Italy and Spain having been hit particularly hard by the first wave in late March and April 2020. The impact of COVID-19 on mortality in our selection of countries seems to be smallest in Germany, although we observe an increase at higher ages near the end of 2020. We refer to Islam et al. (Reference Islam, Shkolnikov, Acosta, Klimkin, Kawachi, Irizarry, Alicandro, Khunti, Yates, Jdanov, White, Lewington and Lacey2021) for a more detailed evaluation of excess mortality in 29 countries, where the baseline mortality
$\hat{D}_{x,\;t,\;w}^i$
in (2) is calculated based on an over-dispersed Poisson regression model accounting for seasonal variation and temporal trends.
For mortality modelling applications such as demographic projections or many actuarial calculations, the development on a yearly scale is of greater importance than weekly fluctuations. One has to note that we only have one complete year of data influenced by COVID-19 at the moment and that it is not clear for how much longer the pandemic will continue to have a strong impact on mortality. However, 2020 marks a notable increase in the age-standardised yearly death rates
$m_{S,\;t}^i$
of some populations such as Belgium, Spain and Poland compared to the previous slightly but steadily decreasing trend. Looking at age-specific rates at higher ages, we have found that the increase of mortality can be even more substantial. One could say that COVID-19 sets back all of these countries in their mortality development by several years. For example, the age-standardised death rates of the Polish male and Spanish female populations in 2020 are at a similar level as in 2008, making them, at least temporarily, lose 12 years of mortality development. The smallest “loss” of development in this sense was suffered by German females, who had higher age-standardised death rates in 2015 than in 2020. As an illustration, the mortality development of these three countries over the last 15 years is depicted in Figure 3.
Yearly German, Polish and Spanish age-standardised death rates
$m_{S,\;t}^i$
as defined in (4) between 2006 and 2020 for females (red, solid), males (green, dash) and the total population (blue, long dash).
We look at age-standardised improvement rates
$I_{S,\;t}^i$
for a comparison in the time dimension. Table 2 lists the 10 worst years in terms of mortality development for the total populations. We see that for all countries 2020 is among the 10 years with the lowest age-standardised improvement rate. Considering only males, 2020 is among the worst 10 years for all the countries as well, whereas it is not among the worst 10 years for Swedish and French females.
Note that the length of the considered time period differs due to limited data availability for some countries. For populations with a long available history of mortality data, we often see that the 1918 pandemic had a larger impact, but apart from that 2020 improvement rates in the considered age range are mostly comparably low to those in years of war or years with other pandemics (such as the 1968–70 Hong Kong flu).
4.2. Mortality shock in the jump-off year
We compare the model trained on real data from 1991 to 2020 (real_20 in Table 1) to the model trained on real data from 1991 to 2019 and a 2020 best estimate (LC_20 in Table 1).
Looking at the model parameters in Figure 4, we find, as predicted by Proposition 1, that the age-dependent parameters
$\alpha_x$
and
$\beta_x$
of both models differ only very slightly. The same holds for the period effects
$\kappa_t$
for
$t < 2020$
. However, for the model which takes the real development of 2020 into account we observe an upward jump in the period effect
$\kappa_{2020}$
for most populations, which might in some cases even lead to a violation of the usual assumption that the period effect increments follow a normal distribution.
Country-specific list of the 10 years with the lowest age-standardised improvement rates
$I_{S,\;t}^i$
as defined in (7). Years are sorted ascendingly by improvement rate so that the worst year is listed first. 2020 is marked in bold.
These differences in calibrated parameters have implications on forecast death rates, in particular because 2020 is the last year in the calibration period so that it is used as the jump-off year for forecasts. This means that it strongly influences the fitted drift of the random walk in (A.2). Unsurprisingly, taking COVID-19 into account leads to higher forecasts. Due to the LC projection methodology, these forecasts drift away from the forecasts by the model which ignores COVID-19 with an increasing forecasting horizon.
We quantify this effect by considering present values of a 30-year annuity for a life aged 65 and a 30-year life insurance contract for a life aged 35 at the beginning of 2021 in Figure 5. Only results for males are shown as the described effects are visible more easily for them, but the picture is qualitatively similar for the female and total populations. Looking at the values based on central forecasts, there are significant differences between countries, with France, Germany and Sweden experiencing the smallest changes in annuity and life insurance policy values. At the other extreme, annuity values for Polish males drop by 9% due to the inclusion of COVID-19, and the value of life insurance policies for Italian males increases by 29%.Footnote
2
Furthermore, we observe that the COVID-19 mortality jump has led to an increase in prediction uncertainty and a corresponding widening of interval forecasts. For the annuity, this is most pronounced for Spanish females with an increase in interval width by a factor of
$2.10$
. For the term assurance, Polish males experience the strongest increase in interval width (factor
$2.35$
).
Country-specific LC model parameters for males, comparing an LC model trained on real data up to 2020 (blue triangle) and an LC model trained on real data up to 2019 and 2020 best estimates (red circle).
Country-specific annuity and life insurance values for males (based on point and interval death rate forecasts), comparing an LC model trained on real data up to 2020 (blue triangle) and an LC model trained on real data up to 2019 and 2020 best estimates (red circle). Discount factor
$v = \frac{1}{1.005}$
.
4.3. Mortality shock before the jump-off year
We compare the model trained on real data from 1992 to 2020 and a best estimate for 2021 (real_20/LC_21 model in Table 1) to the model trained on real data from 1992 to 2019 and 2020/2021 best estimates (LC_20/21 model in Table 1). As mentioned in Section 3.2, best estimates are obtained from an auxiliary LC model calibrated only on data up to 2019 (LC_aux in Table 1).
As expected, these two models produce very similar point forecasts, which usually differ by at most 5%. However, there are material differences between the variance estimates of the associated random walks. The variances are estimated to be up to
$6.7$
times higher for the LC model calibrated on the real 2020 data, which of course results in substantially wider prediction intervals.
This is also evident from Figure 6, where we display annuity and life insurance policy values and their intervals. While the point forecasts are always on a very similar level, the intervals can be much wider when we use the real death rates observed in 2020 instead of the LC central forecast. In the most extreme case (Spanish females), the intervals widen by a factor of
$2.58$
for annuities and
$2.61$
for term assurance. This illustrates that an appropriate way of adjusting for mortality shocks in modelling is highly relevant for applications in which extreme events play a role, for example, in risk management.
Country-specific annuity and life insurance values for males (based on point and interval death rate forecasts), comparing an LC model trained on real data up to 2020 and 2021 best estimates (blue triangle) and an LC model trained on real data up to 2019 and 2020/2021 best estimates (red circle). Discount factor
$v = \frac{1}{1.005}$
.
It is important to note that our observations regarding the widening of prediction intervals do not depend on the exact timing of the observed mortality jump. In particular, the mortality changes induced by COVID-19 will continue to influence uncertainty estimations of the LC model for a long time, unless the model is adjusted in some way. For example, if we were to wait for 10 years and recalibrate the model on data from 2002 to 2031 (assuming a “normal” mortality development between 2021 and 2031), our findings would not change much because the highly volatile years from 2019 to 2021 would still be part of the training data set.
It must further be emphasised that we have assumed that mortality will return to “normal” (LC) levels in 2021 for illustrative purposes in the above calculations. It is possible and, in view of the situation in many countries during the first months of 2021, rather probable that this will not be the case. As an alternative scenario, we have considered the more pessimistic possibility that death rates do not change at all from 2020 to 2021. In this case, the results of comparing an LC model calibrated on these data to a model trained on LC best estimates for 2020 and 2021 are very similar to those obtained in Section 4.2: different point forecasts are produced because the death rates used for calibration differ in the jump-off year, and prediction intervals of the model calibrated on the data including COVID-19 mortality can be substantially wider because the training data contain the mortality shock in 2020.
4.4. Robustness check: using a different calibration method
As detailed in Appendix A, we assume logarithmic death rates to be normally distributed and, therefore, calibrate the LC model by singular value decomposition (SVD). Alternatively, there are also arguments in favour of making a Poisson distribution assumption for the death counts and performing calibration by maximum likelihood estimation (MLE), see Brouhns et al. (Reference Brouhns, Denuit and Vermunt2002). To check the robustness of our results with respect to the calibration method, we have reiterated the above evaluations when calibrating the LC model by MLE instead of SVD.
The model parameters, which are shown in Figure D.1 (Appendix D), are similar to the ones obtained by SVD. For some populations, such as Polish, Italian and Spanish males,
$\kappa_{2020}$
shows a slightly stronger increase as a result of the COVID-19 mortality shock. This also leads to somewhat more pronounced changes in annuity and term assurance values. A comparison of the models real_20 to LC_20 (jump-off year 2020, see Table 1) is displayed in Figure D.2. Annuity values decrease by up to 11% and interval widths increase by a factor of up to
$2.14$
(Polish males), whereas term assurance values increase by up to 44% (Italian males) and interval widths by a factor of up to
$2.90$
(Polish males). With 2021 as the jump-off year and under the assumption of a return to normal (LC best estimate) mortality levels in 2021, the COVID-19 mortality shock has very little influence on point forecasts, see Figure D.3. However, interval widths strongly increase, by a factor of up to
$2.96$
for annuities and up to
$3.08$
for life insurance (Polish males).
Empirically, the two calibration methods lead to qualitatively similar findings, with COVID-19 having an even larger impact under MLE calibration. Our conclusions are robust with respect to the calibration method.
4.5. Robustness check: the Cairns-Blake-Dowd model
To further check the robustness of the results in Sections 4.2 and 4.3, we have repeated our analysis for the CBD model (19). The model parameters
$\kappa_t^{(1)}$
and
$\kappa_t^{(2)}$
are calibrated by linear regression (for every t, regress
$\log m_{x,\;t}$
on
$x-\bar{x}$
), and forecasts are obtained via a two-dimensional random walk with drift (Villegas et al., Reference Villegas, Kaishev and Millossovich2018). Due to the more stringent parametric shape of the CBD model, we have restricted the analysis to the age range
$x_1 = 60, x_2 = 65, \dots, x_7 = 90+$
.
Regarding the model parameters, which are displayed in Figure D.4 (Appendix D), we see an increase in
$\kappa_{2020}^{(1)}$
when using the real 2020 death rates instead of their best estimates for model calibration, which is expected due to the interpretation of this parameter (general mortality level). In a few populations (e.g., Polish males and Swiss males), we also observe a slight increase in
$\kappa_{2020}^{(2)}$
, which indicates that COVID-19 has not led to a perfect parallel shift over the whole age range in these populations, hitting older age groups slightly harder (see Proposition 2). With 2020 as the jump-off year, annuity values drop by up to 13% and interval widths increase by a factor of up to
$1.90$
(Polish males) as a result of the COVID-19 mortality shock, see Figure D.5. With 2021 as the jump-off year and under the assumption of a return to normal (CBD best estimate) mortality levels in 2021, annuity values do not change, but interval widths increase by a factor of up to
$2.76$
(Polish males), see Figure D.6.
Summarising, the results based on the CBD model are qualitatively similar to the ones obtained for the LC model. This illustrates some robustness of our findings to the choice of model.
5. How to Adjust for Mortality Shocks
The findings of the preceding section show that a mortality shock such as the one induced by the COVID-19 pandemic can strongly influence mortality models and the resulting valuation of insurance products. This leads to the question what measures should be taken for jump years like 2020 in the modelling process. There are a lot of different options, from treating these years as completely regular data points to removing COVID-19-related deaths from the data as far as possible, for example, by relying on the log-parallel shift assumption as proposed in Cairns et al. (Reference Cairns, Blake, Kessler and Kessler2020) or by subtracting all excess deaths as proposed in Continuous Mortality Investigation (2020) and in this sense “ignoring” the pandemic in future models and forecasts.
Similarly, outlier analysis can be performed on the fitted period effects
$\left(\hat{\kappa}_t\right)_{t\in\mathcal{T}}$
of the LC model in order to identify extreme mortality changes and remove their influence on the model. Lee & Carter (Reference Lee and Carter1992) themselves follow this idea, but they somewhat arbitrarily only identify the 1918 Spanish flu as an outlier and then apply an intervention model to remove its effect. Li & Chan (Reference Li and Chan2005, Reference Li and Chan2007) propose a more systematic approach applying established techniques from time series outlier analysis to
$\left(\hat{\kappa}_t\right)_{t\in\mathcal{T}}$
. They specify four different types of possible outliers, for example, additive outliers which only affect the series at one point in time, and use t-statistics on the empirical residuals of the assumed time series model to identify them. Finally, they adjust the period effect time series for the detected outliers and recalibrate the time series model.
Outlier adjustment methods restore the (approximate) normality of period effect increments, which means that the random walk model is valid again for the adjusted data. However, as Li & Chan (Reference Li and Chan2007) and Chen & Cox (Reference Chen and Cox2009) note, ignoring mortality jumps might not be applicable, for example, when calculating solvency capital requirements or pricing mortality-linked securities. This is confirmed by our empirical results, which indicate that completely excluding 2020 data from the calibration or applying a similar outlier adjustment method could lead to an underestimation of prediction uncertainty. For instance, the model does in this case not allow for the occurrence of new COVID-19 waves or, generally, further pandemics, even though it is quite probable that they will be a recurring problem in the future (Taubenberger & Morens, Reference Taubenberger and Morens2010; Jonas, Reference Jonas2013; Dodds, Reference Dodds2019; Engel & Ziegler, Reference Engel and Ziegler2020).
On the other hand, treating 2020 as a regular data point by continuing to use the random walk with drift might not be appropriate, either, as it leads to very pessimistic point forecasts and as it is equivalent to a normal distribution assumption for the period effect increments, which may be violated as a result of the COVID-19 mortality shock. However, this normality assumption can be relaxed by a broad range of techniques often used in similar situations in mathematical finance, such as regime switching, extreme value theory or jump processes.
Milidonis et al. (Reference Milidonis, Lin and Cox2011) propose a Markov regime switching model with two hidden states for the LC period effects. They assume that
$\Delta \hat{\kappa}_t \;:\!=\; \hat{\kappa}_t - \hat{\kappa}_{t+1}$
switches between two normal distributions, i.e.,
\begin{equation} \Delta \hat{\kappa}_t \sim \begin{cases} \mathcal{N}\left(\mu, \sigma^2 \right) \mbox{if } \rho_t = 1 \\ \mathcal{N}\left(m, s^2 \right) \mbox{if } \rho_t = 2 \end{cases},\end{equation}
where
$\rho_t$
denotes an underlying, hidden binary Markov chain. Estimates for the transition probabilities
$\mathbb{P}\left(\rho_t = 1|\rho_{t-1} = 2\right)$
and
$\mathbb{P}\left(\rho_t = 2|\rho_{t-1} = 1\right)$
are calibrated by maximum likelihood, along with the remaining parameters
$\mu, m, \sigma, s$
. The model is in principle easily generalised to more than two regimes, and it contains the ordinary LC model as the special case with only one regime. The advantage of regime switching models consists in their precise description of structural changes in the underlying time series, including the possibility for serial correlation in the occurrence of “shock years”. For instance, it is easy to model an extreme event which potentially has an increasing influence on volatility for several years by moving into a high-volatility regime for this duration. The obtained classification of historical mortality observations into two regimes also facilitates interpretation. For example, on 1901–2005 US data, Milidonis et al. (Reference Milidonis, Lin and Cox2011) find a high-volatility regime containing the years around the Spanish flu and a low-volatility regime containing most other years. Drawbacks of this model are its relatively high number of parameters and that it tends to generate rather wide prediction intervals because of the two distinct states which can have very different drift values
$\mu$
and m.
Extreme value theory seeks to model events which occur seldom but can have a large impact if they occur. Chen & Cummins (Reference Chen and Cummins2010) use the peaks-over-threshold approach from extreme value theory to include negative mortality jumps (improvements) in the LC model. They determine a threshold for
$\Delta \hat{\kappa}_t$
above which it is modelled by the usual normality assumption, whereas the tail below the threshold is approximated by a generalised Pareto distribution (GPD). Of course, this mixture distribution can easily be adjusted to focus on upward mortality jumps instead. The peaks-over-threshold approach has deep theoretical foundations, in particular the theorem by Pickands (Reference Pickands1975) and Balkema & de Haan (Reference Balkema and de Haan1974) stating that, under certain conditions, the excess value of a random variable over some threshold converges in distribution to a GPD as the threshold grows. The main difficulty in the calibration of this model lies in identifying a suitable threshold, which ideally requires some domain knowledge regarding the point at which values of
$\Delta\hat{\kappa}_t$
are considered extreme. As a simpler alternative, Chen & Cummins (Reference Chen and Cummins2010) choose it by profile likelihood maximisation.
Another idea is to include a jump process in the time series model for the LC period effects. Jumps can be transitory or permanent (Cox et al., Reference Cox, Lin and Pedersen2010) and their severity can be assumed to follow different distributions (e.g., normal, truncated normal, Pareto, beta). Chen & Cox (Reference Chen and Cox2009) propose an extension of the LC period effect model (A.2) with transitory jumps via
\begin{equation} \hat{\kappa}_{t+1} = \hat{\kappa}_t + \mu + e_{t+1} + N_{t+1} Y_{t+1} - N_t Y_t\end{equation}
where the jump indicator
$N_t$
follows a Bernoulli distribution and the jump severity
$Y_t$
is normally distributed. Deng et al. (Reference Deng, Brockett and MacMinn2012) note that the implicit symmetry assumption of the normal distribution might not be suitable for mortality modelling because upward jumps (mortality deterioration) typically have lower frequency but higher severity than downward jumps (mortality improvement). Inspired by the option pricing model of Kou & Wang (Reference Kou and Wang2004), they propose to account for this by using a double exponential distribution. In our notation, this means that
$\log\left(Y_t\right)$
in (23) is a random variable with density function
\begin{equation} f(y) = p \eta_u e^{-\eta_u y} {\unicode[Times]{x1D7D9}}_{\{y\geq0\}} + q\eta_d e^{\eta_d y} {\unicode[Times]{x1D7D9}}_{\{y<0\}}\end{equation}
where
${\unicode[Times]{x1D7D9}}$
denotes an indicator function, p and q describe the proportion and
$\eta_u,\eta_d$
the severity of upward and downward jumps, respectively. Further generalisations such as age-specific severity patterns (Liu & Li, Reference Liu and Li2015) are possible. Although longer-lasting shocks can occur under the jump model, its conceptual idea is to model transitory shock events which only influence mortality for one time period (1 year). This is suitable for short-term shock events, but might necessitate some adjustment for the COVID-19 mortality shock as it continues to evolve in 2021 and, possibly, in the years to come.
Which of the above options should be chosen depends on the views of the modeler regarding the mentioned advantages and drawbacks and, in particular, on the application. Defining several possible scenarios and comparing the performance of different models under each scenario can aid an informed decision for the most appropriate way of handling extreme mortality events in the modelling process. For this purpose, projections of COVID-19 death counts could be used, which are available from various sources such as the Institute for Health Metrics and Evaluation (2021a). We leave the practical implementation of such an analysis along with more detailed recommendations regarding model choice for future work.
6. Conclusion
We have provided a descriptive analysis of mortality shocks due to COVID-19 in different populations along with a comparison to historical observations. These shocks can have a substantial impact on mortality forecasts and present values of life insurance products. More precisely, our analysis illustrates that a significant impact of a mortality shock on LC point forecasts obtained under the usual random walk assumption is only visible when the shock occurs in the forecast jump-off year. However, we have further observed that prediction uncertainty can increase dramatically whenever a mortality shock appears in the calibration period of the LC model.
Although our quantitative results come with the caveat that official data are still subject to change, for example, due to delayed or erroneous reporting, we are confident that our qualitative conclusions are valid. Nevertheless, it would be interesting to repeat the analysis as more recent data become available. Looking at the first months of 2021, the hypothetical scenario of a quick return to normal mortality levels which we consider in Section 4.3 is unrealistic in many countries. This means that our main conclusions and recommendations for 2020 might turn out to be applicable to 2021 as well.
Many questions still remain open, in particular regarding the future development of COVID-19-related mortality:
-
• SARS-CoV-2 is generally expected to become endemic (Phillips, Reference Phillips2021). However, the recent uptake in vaccinations is believed to substantially reduce COVID-19-related mortality (Ng & Reid, Reference Ng and Reid2021; Andrews et al., Reference Andrews, Keeling, Stowe, Ismael, Coughlan, Allen, Ramsay and Lopez Bernal2021; Moore et al., Reference Moore, Hill, Tildesley, Dyson and Keeling2021). Given these and other influence factors such as ongoing non-pharmaceutical interventions, which pattern of COVID-19 deaths will we observe in the future?
-
• Will COVID-19 cause new cohort effects, for example, selection effects because primarily the healthier individuals might survive whereas the more frail might perish (Cairns et al., Reference Cairns, Blake, Kessler and Kessler2020)? It would be interesting to repeat our analysis with a model including cohort effects to see whether COVID-19 leads to changes in the cohort parameters and, if yes, which cohorts are influenced the most.
-
• Will there be amplifying or weakening interaction effects between COVID-19 and Influenza?
While these questions are impossible to answer definitively, they give rise to several possible scenarios that could be investigated in future research.
Finally, the most important question for mortality modelers right now is of course how 2020 – and possibly also 2021 – mortality data should be treated in the modelling process. We have presented several approaches in Section 5. Further research should focus on evaluating their suitability for current mortality data and deriving specific, problem-adjusted recommendations for practitioners. Ideally, a model adjustment technique can be identified which produces satisfactory results over a broad range of scenarios.
Data Availability Statement
All data used in the paper are publicly available at https://mortality.org. Code for reproducing the results of the paper is publicly available at https://github.com/schnuerch/COVID-Impact.
Acknowledgments
We thank two anonymous reviewers for their thoughtful suggestions, which have improved the manuscript. Furthermore, we thank Renate Ernst for her helpful advice. S. Schnürch is grateful for the financial support from the Fraunhofer Institute for Industrial Mathematics ITWM. T. Kleinow acknowledges financial support from the Actuarial Research Centre of the Institute and Faculty of Actuaries through the research programme on “Modelling, Measurement and Management of Longevity and Morbidity Risk”.
Appendix A. Lee-Carter Model Fitting and Forecasting
For a fixed population i, Lee & Carter (Reference Lee and Carter1992) model logarithmic death rates as
\begin{equation} \log m_{x,\;t}^i = \alpha_x^i + \beta_x^i \kappa_t^i + \varepsilon_{x,\;t}^i\end{equation}
where
$\alpha_x^i$
is an age-specific base mortality level,
$\kappa_t^i$
is the period effect which describes mortality development over time and
$\beta_x^i$
is the age effect which allows changes in the period effect to influence different ages with varying intensity. In our numerical studies in Section 4, we consider data in 5-year age groups, which means, for example, that
$\beta_{35}$
denotes the age effect for the age group
$35-39$
. The error terms
$\varepsilon_{x,\;t}^i$
are usually assumed to follow homoskedastic normal distributions. The model can then be fit via singular value decomposition. In the following, we suppress the dependence on i in our notation.
It is easily seen that the parameters of the LC model are not identifiable without imposing some constraints. Several types of identifiability constraints have been used by different authors (Hunt & Blake, Reference Hunt and Blake2020). Here, we will use (14), i.e.,
$\sum_{x=x_1}^{x_A} \beta_x = 1$
and
$\kappa_{t_1} = 0$
. It is always possible to transform parameters obeying (14) to parameters obeying any other set of identifiability constraints, and vice versa. The choice of identifiability constraints has by definition no influence on the fitted mortality rates. It can have an influence on point and distribution forecasts, which can be avoided by using a location and scale invariant forecasting technique such as an ARIMA model with intercept, see Nielsen & Nielsen (Reference Nielsen and Nielsen2014), Hunt & Blake (Reference Hunt and Blake2020) for details.
To obtain mortality forecasts, Lee & Carter (Reference Lee and Carter1992) propose to fit such an ARIMA model to the calibrated time series of period effects
$\left(\hat{\kappa}_t\right)_{t\in\mathcal{T}}$
. The standard choice is a random walk with drift, i.e.,
\begin{equation} \hat{\kappa}_{t+1} = \mu + \hat{\kappa}_t + e_{t+1} \text{ with } \mu\in\mathbb{R} \text{ and } e_{t+1} \sim \mathcal{N}\left(0, \sigma^2\right) \text{ i.i.d.}\end{equation}
Point forecasts for year
$t_Y+h$
are then obtained by inserting the expectation of
$\hat{\kappa}_{t_Y+h}$
in the LC model, i.e.,
\begin{equation} \hat{m}_{x,\;t_Y+h} \;:\!=\; \exp\left(\hat{\alpha}_x + \hat{\beta}_x \left(\hat{\kappa}_{t_Y} + h\cdot \hat{\mu}\right)\right)\end{equation}
with the estimator
\begin{equation} \hat{\mu} \;:\!=\; \frac{\hat{\kappa}_{t_Y} - \hat{\kappa}_{t_1}}{Y-1}\end{equation}
for the drift.
An important advantage of stochastic mortality models is the possibility to obtain estimates of prediction uncertainty in addition to central forecasts, which is, for example, vital in risk management applications. We calculate h-step prediction intervals at level
$a\in(0,1)$
by
\begin{equation} \hat{m}_{x,\;t_Y+h}^{\text{lower}|\text{upper}} \;:\!=\; \exp\left(\hat{\alpha}_x + \hat{\beta}_x\left(\hat{\kappa}_{t_Y} + h\cdot\hat{\mu} \pm \sqrt{h} \cdot\hat{\sigma} \cdot \Phi^{-1}\left(\frac{1 + a}{2}\right) \right)\right)\end{equation}
where
$\Phi^{-1}$
denotes the quantile function of the standard normal distribution and
\begin{equation} \hat{\sigma}^2 \;:\!=\; \frac{1}{Y-2}\sum\limits_{t=t_2}^{t_Y} \left(\hat{\kappa}_t - \hat{\kappa}_{t-1}\right)^2\end{equation}
is the estimator for
$\sigma^2$
. For our numerical results in Section 4, we use
$a=0.95$
, i.e., we calculate 95% prediction intervals.
Note that we do not include parameter uncertainty here to allow for better comparability to other models. With increasing forecasting horizon, the influence of parameter uncertainty on total prediction uncertainty grows and so, depending on the application, it should be taken into account as well. We refer to Kleinow & Richards (Reference Kleinow and Richards2016) for a general treatment of parameter uncertainty estimation in period effect ARIMA models.
Appendix B. Recalibration Consistency of the LC Model
The following proposition shows that the LC model is consistent under recalibration.
Proposition 3 (Recalibration consistency). Consider an LC model calibrated to
$m_{x,\;t}$
where
$x\in\mathcal{X}$
and
$t\in\{t_1,\dots,t_{Y-1}\}$
, and denote its one-step forecasts by
$\hat{m}_{x,\;t_Y}$
.
Consider another LC model calibrated to
$m_{x,\;t}^{\text{new}}$
,
$x\in\mathcal{X}$
and
$t\in\mathcal{T} = \{t_1,\dots,t_Y\}$
, where
$m_{x,\;t}^{\text{new}} \;:\!=\; m_{x,\;t}$
for
$t\in\{t_1,\dots,t_{Y-1}\}$
and
$m_{x,\;t_Y}^{\text{new}} \;:\!=\; \hat{m}_{x,\;t_Y}$
.
These two models yield the same fitted values
$\alpha_x + \beta_x \kappa_t$
for all
$x\in\mathcal{X}$
and
$t\in\mathcal{T}$
and so can only differ in the choice of identifiability constraints. In particular, using a scale and location invariant forecasting method in the sense of Hunt & Blake (Reference Hunt and Blake2020), their forecasts beyond
$t_Y$
also exactly coincide.
Proof. This follows from a more general result in the setting of learning an estimator by empirical risk minimisation. Given some training data set
$(o_1, y_1), \dots, (o_n, y_n)$
with observations
$o_1,\dots,o_n$
and ground truths
$y_1,\dots,y_n$
, the empirical risk of an estimator h is defined as
\begin{equation} \hat{R}_n(h)\;:\!=\; \frac{1}{n} \sum\limits_{j=1}^n \ell\left(y_j, h(o_j)\right) \end{equation}
for some real-valued loss function
$\ell$
. Denoting a minimiser of (B.1) by
$h^{\text{min}}$
and assuming that another data point
$(o_{n+1}, h^{\text{min}}(o_{n+1}))$
is added to the training data set, the adjusted empirical risk
\begin{equation} \hat{R}_{n+1}(h) = \frac{1}{n+1} \left(\sum\limits_{j=1}^n \ell\left(y_j, h(o_j)\right) + \ell\left(h^{\text{min}}(o_{n+1}), h(o_{n+1}) \right)\right) \end{equation}
is clearly still minimised by
$h^{\text{min}}$
, i.e., the model does not change (up to identifiability issues).
It is easy to see that calibrating an LC model amounts to minimising (B.1) with a suitable loss function, for example, quadratic loss under formulation (A.1), and the claim follows.
While the fit stays unchanged, the uncertainty estimates (A.5) resulting from the two models described in Proposition 3 would differ.
Appendix C. Proofs of Results from Section 3.3
Proof of Proposition 1. We calculate
\begin{equation*} \begin{aligned} \alpha_x &= \log m_{x,\;t_1} - \beta_x \kappa_{t_1} - \varepsilon_{x,\;t_1} = \log m_{x,\;t_1}^B - \varepsilon_{x,\;t_1} = \alpha_x^B - \Delta \varepsilon_{x,\;t_1} \end{aligned} \end{equation*}
which shows (15),
\begin{equation*} \begin{aligned} \kappa_t &= \sum_{x=x_1}^{x_A} \beta_x \kappa_t = \sum_{x=x_1}^{x_A} \left(\log m_{x,\;t} - \alpha_x - \varepsilon_{x,\;t}\right) \\ &= \sum_{x=x_1}^{x_A} \left(\log m_{x,\;t}^B + \log\left(1+c\right){\unicode[Times]{x1D7D9}}_{t=t_Y} - \alpha_x^B + \Delta\varepsilon_{x,\;t_1} - \varepsilon_{x,\;t}\right) \\ &= \kappa_t^B + A\log\left(1+c\right) {\unicode[Times]{x1D7D9}}_{t=t_Y} + \sum_{x=x_1}^{x_A} \left(\Delta \varepsilon_{x,\;t_1} - \Delta \varepsilon_{x,\;t}\right) \end{aligned} \end{equation*}
which shows (17), and
\begin{equation*} \begin{aligned} \beta_x &= \frac{\log m_{x,\;t^*} - \alpha_x - \varepsilon_{x,\;t^*}}{\kappa_{t^*}} \\ &= \frac{\log m_{x,\;t^*}^B - \alpha_x^B + \Delta\varepsilon_{x,\;t_1} - \varepsilon_{x,\;t^*}}{\kappa_{t^*}^B + \sum_{x=x_1}^{x_A} \left(\Delta \varepsilon_{x,\;t_1} - \Delta \varepsilon_{x,\;t^*}\right)} \\ &= \frac{\beta_x^B\kappa_{t^*}^B + \Delta\varepsilon_{x,\;t_1} - \Delta\varepsilon_{x,\;t^*}}{\kappa_{t^*}^B + \sum_{x=x_1}^{x_A} \left(\Delta \varepsilon_{x,\;t_1} - \Delta \varepsilon_{x,\;t^*}\right)} \approx \frac{\beta_x^B\kappa_{t^*}^B}{\kappa_{t^*}^B} \end{aligned} \end{equation*}
where
$t^*\not\in\{t_1, t_Y\}$
is chosen such that
$\kappa_{t^*} \neq 0$
and
$\kappa_{t^*}^B \neq 0$
, which shows (16).
Proof of Corollary 1. The assumption that the period effect follows a random walk with drift means that Equation (A.2) holds for
$\hat{\kappa}_{t+1}$
and, suitably adjusting the notation, for
$\hat{\kappa}_{t+1}^B$
. The maximum likelihood estimator for the drift
$\mu$
is
\begin{equation} \hat{\mu} \;:\!=\; \frac{\hat{\kappa}_{t_Y} - \hat{\kappa}_{t_1}}{Y-1} \approx \frac{\hat{\kappa}_{t_Y}^B + A\log\left(1+c\right) - \hat{\kappa}_{t_1}^B}{Y-1} = \hat{\mu}^B + \frac{A}{Y-1} \log\left(1+c\right) \end{equation}
where we have used equation (17) from Proposition 1. Combined with equations (15) and (16), this implies for the central forecast as defined in (A.3)
\begin{equation*} \begin{aligned} \hat{m}_{x,\;t_Y+h} &= \exp\left(\hat{\alpha}_x + \hat{\beta}_x \left(\hat{\kappa}_{t_Y} + h\hat{\mu}\right)\right) \\ &\approx \exp\left(\hat{\alpha}_x^B + \hat{\beta}_x^B \left(\hat{\kappa}_{t_Y}^B + A\log\left(1+c\right) + h\left(\hat{\mu}^B + \frac{A}{Y-1} \log\left(1+c\right)\right)\right)\right) \\ &= \hat{m}_{x,\;t_Y+h}^B \cdot \exp\left(\hat{\beta}_{x}^B \cdot A \log\left(1+c\right) \left(1 + \frac{h}{Y-1}\right) \right) \end{aligned} \end{equation*}
Proof of Proposition 2. We calculate
\begin{equation*} \begin{aligned} \kappa_t^{(1)} &= \kappa_t^{(1)} + (\bar{x} - \bar{x})\kappa_t^{(2)} \\ &= \log m_{\bar{x}\;,t} - \varepsilon_{\bar{x},\;t} \\ &= \log m_{\bar{x},\;t}^B + \log\left(1+c\right){\unicode[Times]{x1D7D9}}_{t=t_Y} - \varepsilon_{\bar{x},\;t} \\ &= \kappa_t^{(1),B} + \log\left(1+c\right){\unicode[Times]{x1D7D9}}_{t=t_Y} - \Delta \varepsilon_{\bar{x},\;t} \end{aligned} \end{equation*}
and
\begin{equation*} \begin{aligned} \kappa_t^{(2)} &= \frac{\log m_{x_1,\;t} - \varepsilon_{x_1,\;t} - \log m_{\bar{x},\;t} + \varepsilon_{\bar{x},\;t}}{x_1 - \bar{x}} \\ &= \frac{\log m_{x_1,\;t}^B - \varepsilon_{x_1,\;t} - \log m_{\bar{x},\;t}^B + \varepsilon_{\bar{x},\;t}}{x_1 - \bar{x}} \\ &= \kappa_t^{(2),B} + \frac{\Delta\varepsilon_{\bar{x},\;t} - \Delta\varepsilon_{x_1,\;t}}{x_1-\bar{x}} \end{aligned} \end{equation*}
Appendix D. Additional Figures
Country-specific LC model parameters for males, comparing an LC model trained on real data up to 2020 (blue triangle) and an LC model trained on real data up to 2019 and 2020 best estimates (red circle), calibration method: Poisson maximum likelihood estimation.
Country-specific annuity and life insurance values for males (based on point and interval death rate forecasts), comparing an LC model trained on real data up to 2020 (blue triangle) and an LC model trained on real data up to 2019 and 2020 best estimates (red circle), calibration method: Poisson maximum likelihood estimation. Discount factor
$v = \frac{1}{1.005}$
.
Country-specific annuity and life insurance values for males (based on point and interval death rate forecasts), comparing an LC model trained on real data up to 2020 and 2021 best estimates (blue triangle) and an LC model trained on real data up to 2019 and 2020/2021 best estimates (red circle), calibration method: Poisson maximum likelihood estimation. Discount factor
$v = \frac{1}{1.005}$
.
Country-specific CBD model parameters for males, comparing a CBD model trained on real data up to 2020 (blue triangle) and a CBD model trained on real data up to 2019 and 2020 best estimates (red circle).
Country-specific annuity values
$a_{65:\overline{30}\kern-0.5pt\raise.5pt\hbox{${\scriptstyle{\mid}}$}}(2021)$
for males (based on point and interval death rate forecasts), comparing a CBD model trained on real data up to 2020 (blue triangle) and a CBD model trained on real data up to 2019 and 2020 best estimates (red circle). Discount factor
$v = \frac{1}{1.005}$
.
Country-specific annuity values
$a_{65:\overline{30}\kern-0.5pt\raise.5pt\hbox{${\scriptstyle{\mid}}$}}(2022)$
for males (based on point and interval death rate forecasts), comparing a CBD model trained on real data up to 2020 and 2021 best estimates (blue triangle) and a CBD model trained on real data up to 2019 and 2020/2021 best estimates (red circle). Discount factor
$v = \frac{1}{1.005}$
.
Open access
