# Mortality Forecasting using Factor Models: Time-varying or Time-invariant
Factor Loadings?

Lingyu He  Hunan University  Fei Huang Correspondence to: Dr. Fei Huang,
School of Risk and Actuarial Studies, UNSW Business School, UNSW Sydney, NSW
2052, Australia. Email<EMAIL_ADDRESS>The University of New South Wales
Jianjie Shi  Monash University  Yanrong Yang  The Australian National
University

###### Abstract

Many existing mortality models follow the framework of classical factor
models, such as the Lee-Carter model and its variants. Latent common factors
in factor models are defined as time-related mortality indices (such as
$\kappa_{t}$ in the Lee-Carter model). Factor loadings, which capture the
linear relationship between age variables and latent common factors (such as
$\beta_{x}$ in the Lee-Carter model), are assumed to be time-invariant in the
classical framework. This assumption is usually too restrictive in reality as
mortality datasets typically span a long period of time. Driving forces such
as medical improvement of certain diseases, environmental changes and
technological progress may significantly influence the relationship of
different variables. In this paper, we first develop a factor model with time-
varying factor loadings (time-varying factor model) as an extension of the
classical factor model for mortality modelling. Two forecasting methods to
extrapolate the factor loadings, the local regression method and the naive
method, are proposed for the time-varying factor model. From the empirical
data analysis, we find that the new model can capture the empirical feature of
time-varying factor loadings and improve mortality forecasting over different
horizons and countries. Further, we propose a novel approach based on change
point analysis to estimate the optimal ‘boundary’ between short-term and long-
term forecasting, which is favoured by the local linear regression and naive
method, respectively. Additionally, simulation studies are provided to show
the performance of the time-varying factor model under various scenarios.

JEL Code  G22 - Insurance; Insurance Companies; Actuarial Studies

Keywords  Lee-Carter Model; Long-term Forecasting; Optimal ‘Boundary’
Estimation; Short-term Forecasting; Time-varying Factor Model.

##  1 Introduction

Mortality forecasting is an important topic in various areas, such as
demography, actuarial science and government policymaking. Most age-specific
mortality data are high-dimensional time series. The factor model approach is
one of the most popular methods to model high-dimensional time series,
representing the data matrix by a few latent common factors. Common factors
describe common information shared by cross-sections, while factor loadings
reflect the linear relationship between the original variables and the common
factors. There is a large literature discussing factor models, including but
not limited to Anderson (1963), Pena and Box (1987), Stock and Watson (2002),
Bai and Ng (2002), Bai (2009), Lam and Yao (2012) and Chang et al. (2018).

Many existing stochastic mortality models use the factor model approach. As an
application of the classical factor model (with time-invariant factor
loadings), Lee and Carter (1992) (Lee-Carter Model) is one of the most
prominent methods for mortality forecasting, which is employed by the US
Bureau of the Census as the benchmark model to predict long-run life
expectancy (Hollmann et al., 1999). The common factor extracted by the Lee-
Carter model is defined as Mortality Index, and the factor loadings capture
the relationship between the age variables and the mortality index. Since
there is only one factor in the Lee-Carter model, Booth et al. (2002), Renshaw
and Haberman (2003) and Yang et al. (2010) extended the Lee-Carter framework
to incorporate more common latent factors for mortality modelling in different
countries. Li and Chan (2005) proposed an outlier-adjusted model to deal with
possible outliers in the mortality index by combining the Lee-Carter model
with time series outlier analysis. Additionally, Booth et al. (2006) compared
the Lee-Carter model with four other variants by applying them to mortality
data of multiple populations. Tuljapurkar et al. (2000) examined mortality
rates over five decades for the G7 countries using the Lee-Carter model.
Lundström and Qvist (2004) and Booth et al. (2004) applied the Lee-Carter
model to mortality data of Sweden and Australia, respectively. A summary of
the variants of the Lee-Carter model is discussed in Booth and Tickle (2008).

In the existing literature of mortality factor models, factor loadings, which
capture the relationship between age variables and latent common factors, are
usually assumed to be invariant over time (we call factor models with time-
invariant factor loadings ‘classical factor model’). For example, in Lee-
Carter model, there is only one factor and the time-invariant factor loading
represents the age-related sensitivity to the mortality improvement (we call
classical factor model with only one factor ‘Lee-Carter model’ throughout this
paper). However, since mortality datasets typically span a long period of
time, it is restrictive to assume that the factor loadings are time-invariant.
Driving forces such as medical improvement of certain diseases, environmental
changes, and technological progress may influence the relationship of
different variables significantly. Booth et al. (2002) studied the violation
of the invariance assumption in the mortality data of Australia and suggested
to find an optimal fitting period during which the factor loadings were
invariant to improve the fit of the Lee-Carter model. Their approach, however,
needs to manually select the fitting period and hence loses the information of
early years. In recent years, there is a rich literature on time-varying
factor models to capture the dynamics and structural changes in factor
loadings for macroeconomic variables modelling, for example, see Breitung and
Eickmeier (2011) and Chen et al. (2014). However, there has been no literature
on mortality modelling which allows factor loadings to change smoothly over
time, to the best of our knowledge. Li and O’Hare (2017) and Li et al. (2015)
used semi-parametric approaches to extend the CBD models (Cairns et al., 2009)
by allowing for time-varying coefficients, which can free model assumptions
and show superior short-term forecasting performance. However, CBD models are
only suitable for old-age mortality modelling, and the factors (regressors)
are observable. Unfortunately, for Lee-Carter model and many of its variants
the factors are unobserved, which makes it difficult to model and estimate. To
fill those gaps, we introduce a factor model with time-varying factor loadings
as an extension of the classical factor model based on Su and Wang (2017).
This new model can be used for mortality modelling and forecasting by
developing corresponding estimation and forecasting methods.

As the time-varying factor model allows for time-varying factor loadings, it
provides more flexibility in model fitting, which, however, also poses
challenges in model forecasting. Besides forecasting the common factors,
factor loadings also need to be extrapolated into the future. In this paper,
we provide two forecasting methods of the factor loadings, one uses the local
linear regression to roll over the time-varying factor loadings into the
future; while the other one inherits the value of the factor loading from the
last time period and remains invariant in the future. These two forecasting
methods are called the local linear regression and the naive method,
respectively. Their details are described in Section 2. Empirical results
using the mortality data from different populations show that the time-varying
factor model provides more accurate out-of-sample forecasting results than the
Lee-Carter model.

The existing literature suggests that different forecasting horizons may
favour different models. For example, Bell (1997) found that a simple random
walk with drift model for age-specific mortality rates yields the most
accurate 1-step-ahead forecast compared with the other six methods on the US
data. Hyndman and Ullah (2007) introduced a method which outperformed the
method proposed by Lee and Miller (2001) in the long-term forecasting.
Specifically, we have found in the literature that semi-parametric or non-
parametric methods can be more suitable for short-term forecasting. For
example, the semi-parametric model developed in Li et al. (2015) can produce
superior 5-year-ahead forecasting results. CMI (2009) employed the P-splines
model (Currie et al. (2004)) for short-term forecasting to generate the
initial rates of mortality improvement. Our empirical applications in Section
5.3 also suggest that the time-varying model based on local regression (non-
parametric forecasting) is better for short-term forecasting, while the time-
varying model based on naive method (parametric forecasting) is better for
long-term forecasting. Then where is the optimal ‘boundary’ between short-term
(based on the local regression method) and long-term (based on the naive
method) forecasting? We propose a novel approach based on change point
analysis (Bai, 2010) to estimate the optimal ‘boundary’ and apply it to
mortality data of multiple countries. Additionally, we conduct simulation
studies to show the performance of the time-varying factor model under
different scenarios and investigate under which conditions it preforms better
than the classical factor model.

The rest of the paper is organized as follows. Section 2 introduces the time-
varying factor model and its estimation approach. The forecasting methods
based on the time-varying factor model are also discussed in detail. Section 3
discusses the relative advantages of the local regression method and the naive
method in the short-term and long-term forecasting, respectively. We then
propose an approach based on change point analysis to estimate the ‘boundary’
between short-term and long-term forecasting, which is favoured by the local
regression method and the naive method, respectively. Section 4 introduces the
datasets and empirical evidence of time-varying factor loadings. Section 5
applies the proposed methods to age-specific mortality data of multiple
countries and shows the advantages of the proposed methods. Section 6 conducts
simulation studies to investigate the performance of the time-varying factor
model under different scenarios. Section 7 concludes the paper. Appendix A
provides the gender-specific empirical results using the time-varying factor
model. Appendix B presents the estimations of the optimal boundaries for
multiple countries with a variety of forecasting horizons. Appendix C displays
estimation results of the time-varying model with multiple factors.

##  2 Time-varying Factor Model

Let $m_{x,t}$ denote the central death rate for age $x$ in year $t$, where
$x=1,2,\dots,N$ and $t=1,2,\dots,T$. Thus,
$\\{m_{x,t}\\}_{x=1,2,\dots,N,t=1,2,\dots,T}$ is an $N$-dimensional time
series with $T$ observations. Since mortality rates are always positive
numbers, we use the log transformation to map the central death rates from
$\mathbb{R}^{+}$ space to $\mathbb{R}$ space for modelling purposes. Assume
$a_{x}$ is the age-specific constant, which is the averages over time of the
$\ln(m_{x,t})$. Then $\ln(m_{x,t})-a_{x}$ can be modelled using the classical
factor model as follows:

$\displaystyle\ln(m_{x,t})=a_{x}+{\bm{b}}_{x}^{\top}{\bm{k}}_{t}+\varepsilon_{x,t},$
(2.1)

where ${\bm{k}}_{t}$ is an $R\times 1$ vector of common factors;
${\bm{b}}_{x}$ is an $R\times 1$ vector of factor loadings, capturing the
impact of each common factor on age $x$ (i.e. the age-related sensitivity to
the mortality improvement); and $\varepsilon_{x,t}$ is the idiosyncratic error
of $\ln(m_{x,t})$, which represents the component unexplained by the common
factor. Here, ${\bm{k}}_{t}$, ${\bm{b}}_{x}$ and $\varepsilon_{x,t}$ are all
unobservable components. Specifically, when $R=1$, the classical factor model
is equivalent to the Lee-Carter model. The single factor $k_{t}$ is defined as
the mortality index in the Lee-Carter model, and consequently the factor
loading $b_{x}$ represents the impact of the mortality index on the death rate
of age $x$.

The classical factor model, however, is too restrictive when used to analyse
the mortality data. It assumes that for each age the factor loadings are
invariant over time. Statisticians and economists have noticed that the
relationship between many economic variables and common factors is not time-
invariant. Our empirical analysis using mortality data in Section 5 also
suggests time-varying factor loadings. Therefore, we develop a factor model to
allow for factor loadings changing smoothly overtime.

We introduce the time-varying factor model based on the work of Su and Wang
(2017), where factor loadings are modelled as non-random functions of time. Su
and Wang (2017) provided a localized PCA method to consistently estimate the
factors and time-varying factor loadings. Compared with Park et al. (2009),
the time-varying factor model proposed by Su and Wang (2017) can capture more
types of structural changes in factor loadings, including both continuous
changes and abrupt structural breaks. Assume $\ln(m_{x,t})-a_{x}$ follows the
time-varying factor model with $R$ unobservable common factors:

$\displaystyle\ln(m_{x,t})=a_{x}+{\bm{b}}_{x,t}^{\top}{\bm{k}}_{t}+\varepsilon_{x,t},$
(2.2)

where notations above are the same as the classical factor model, except for
the factor loadings. Here, each component of the factor loading
${\bm{b}}_{x,t}$ is assumed to be a deterministic function of $t/T$:
${\bm{b}}_{x,t}={\bm{b}}_{x}(t/T)$, where each component of
${\bm{b}}_{x}(\cdot)$ is an unknown piece-wise smooth function of $t/T$. The
time-varying factor model can be seen as a generalization of the classical
factor model. If ${\bm{b}}_{x}(\cdot)$ is invariant over time, which is a
special case of the piece-wise smooth function, the time-varying factor model
will degenerate to the classical factor model. Generally speaking, the
assumption that factor loadings are time-invariant is too restrictive to hold
in most settings. However, the time-varying factor model can relax this
assumption by allowing for both continuous structural changes and abrupt
changes in factor loadings, which can also benefit the mortality forecasting.

###  2.1 Identification Problem

Similar to the classical factor model, there exists an identification problem
in the time-varying factor model. At each time point $t$, and for any $R\times
R$ invertible matrix ${\bm{H}}_{t}$, we have
${\bm{b}}_{x,t}^{\top}{\bm{k}}_{t}=\left({\bm{H}}_{t}^{-1}{\bm{b}}_{x,t}\right)^{\top}\left({\bm{H}}_{t}^{\top}{\bm{k}}_{t}\right)$.
Since an $R\times R$ invertible matrix has $R^{2}$ free elements, $R^{2}$
restrictions are needed in parameter estimations so that ${\bm{b}}_{x,t}$ and
${\bm{k}}_{t}$ can be identified separately. Define
${\bm{B}}_{t}=\left({\bm{b}}_{1,t},{\bm{b}}_{2,t},\dots,{\bm{b}}_{N,t}\right)^{\top}$
and
${\bm{K}}=\left({\bm{k}}_{1},{\bm{k}}_{2},\dots,{\bm{k}}_{T}\right)^{\top}$.
Then the two sets of restrictions to solve the issue of identification are
${\bm{K}}^{\top}{\bm{K}}/T=\mathbb{I}_{R}$ and
${\bm{B}}_{t}^{\top}{\bm{B}}_{t}=\text{a diagonal matrix}$, where
$\mathbb{I}_{R}$ is an $R\times R$ identity matrix. The first normalization
condition imposes $R\times(R+1)/2$ restrictions on the parameters, and the
remaining $R\times(R-1)/2$ restrictions are obtained by requiring the second
constraint. These restrictions can uniquely determine the factors ${\bm{K}}$
and the factor loadings ${\bm{B}}_{t}$ (only up to a sign change, i.e.,
$-{\bm{K}}$ and $-{\bm{B}}_{t}$ also satisfy the two sets of restrictions).
When $R=1$, only one restriction is needed to identify parameters. We choose
to use the same normalization condition as Lee and Carter (1992), that is, we
normalize the ${\bm{b}}_{x,t}$ to sum to unity for each $t$. In this way, we
can directly compare the results of our new method with that of the Lee-Carter
model.

###  2.2 Estimation Method

The estimation method for the time-varying factor model is proposed by Su and
Wang (2017). Let $r\in\\{1,\dots,T\\}$ be a fixed year. Since we have assumed
that each component of ${\bm{b}}_{x,t}:[0,1]\rightarrow\mathbb{R}$ is a piece-
wise smooth function, we have:

$\displaystyle{\bm{b}}_{x,t}={\bm{b}}_{x}(\frac{t}{T})\approx{\bm{b}}_{x}(\frac{r}{T})={\bm{b}}_{x,r},\
\text{when}\ \frac{t}{T}\approx\frac{r}{T}.$

Thus, the mortality rate $\ln(m_{x,t})$ can be approximated by:

$\displaystyle\ln(m_{x,t})\approx
a_{x}+{\bm{b}}_{x,r}^{\top}{\bm{k}}_{t}+\varepsilon_{x,t},\ \text{when}\
\frac{t}{T}\approx\frac{r}{T}.$

In order to estimate the factors and time-varying factor loadings, we consider
the following local weighted least squares problem:

$\displaystyle\min_{\\{{\bm{b}}_{x,r}\\}_{x=1}^{N},\\{{\bm{k}}_{t}\\}_{t=1}^{T}}(NT)^{-1}\sum_{x=1}^{N}\sum_{t=1}^{T}\left(\ln(m_{x,t})-a_{x}-{\bm{b}}_{x,r}^{\top}{\bm{k}}_{t}\right)^{2}K_{h}\left(\frac{t-r}{T}\right),$
(2.3)

subject to the identification constraints as discussed in Section 2.1. In the
objective function in Equation (2.3), $K_{h}(x)=h^{-1}K(x/h)$, where
$K\left(\cdot\right)$ is a kernel function and $h$ is a smoothing parameter
called “bandwidth”. We will show that the optimization problem of Equation
(2.3) can be solved using the same estimation method for the classical factor
model.

We have known that the mortality rates can be approximated by
$\ln(m_{x,t})-a_{x}\approx{\bm{b}}_{x,r}^{\top}{\bm{k}}_{t}+\varepsilon_{x,t}$
when $\frac{t}{T}\approx\frac{r}{T}$. Multiplying both sides of the equation
by

$\displaystyle
K_{h,tr}^{1/2}\coloneqq\left(K_{h}\left(\frac{t-r}{T}\right)\right)^{1/2}=\left(\frac{1}{h}K\left(\frac{t-r}{Th}\right)\right)^{1/2},$

we obtain a transformed model as:

$\displaystyle K_{h,tr}^{1/2}(\ln(m_{x,t})-a_{x})\approx
K_{h,tr}^{1/2}{\bm{b}}_{x,r}^{\top}{\bm{k}}_{t}+K_{h,tr}^{1/2}\varepsilon_{x,t},\
\text{when}\ \frac{t}{T}\approx\frac{r}{T}.$

Then we can define matrices

$\displaystyle{\bm{M}}^{(r)}=\left({\bm{M}}^{(r)}_{1},\dots,{\bm{M}}^{(r)}_{N}\right),$
$\displaystyle\bm{\varepsilon}^{(r)}=\left(\bm{\varepsilon}_{1}^{(r)},\dots,\bm{\varepsilon}_{N}^{(r)}\right),$

and

$\displaystyle{\bm{K}}^{(r)}=\left(K_{h,1r}^{1/2}{\bm{k}}_{1},\dots,K_{h,Tr}^{1/2}{\bm{k}}_{T}\right)^{\top}.$

where
${\bm{M}}^{(r)}_{x}=\left(K_{h,1r}^{1/2}(\ln(m_{x,1})-a_{x}),\dots,K_{h,Tr}^{1/2}(\ln(m_{x,T})-a_{x})\right)^{\top}$
and
$\bm{\varepsilon}_{x}^{(r)}=\left(K_{h,1r}^{1/2}\varepsilon_{x,1},\dots,K_{h,Tr}^{1/2}\varepsilon_{x,T}\right)^{\top}$
with $x=1,2,\ldots,N$. Therefore, the transformed model can be written in
matrix form as follows:

$\displaystyle{\bm{M}}^{(r)}\approx{\bm{K}}^{(r)}{\bm{B}}_{r}^{\top}+\bm{\varepsilon}^{(r)},$

and the optimization problem above can also be written in matrix notation as:

$\displaystyle\min_{{\bm{K}}^{(r)},{\bm{B}}_{r}}\operatorname{Tr}{\left(\left({\bm{M}}^{(r)}-{\bm{K}}^{(r)}{\bm{B}}_{r}^{\top}\right)\left({\bm{M}}^{(r)}-{\bm{K}}^{(r)}{\bm{B}}_{r}^{\top}\right)^{\top}\right)}$
$\displaystyle s.t\ {\bm{K}}^{(r)\top}{\bm{K}}^{(r)}/T=\mathbb{I}_{R}\
\text{and}\ {\bm{B}}_{r}^{\top}{\bm{B}}_{r}=\text{a diagonal matrix}.$

Concentrating out
${\bm{B}}_{r}={\bm{M}}^{(r)\top}{\bm{K}}^{(r)}\left({\bm{K}}^{(r)\top}{\bm{K}}^{(r)}\right)^{-1}$
(which is ${\bm{M}}^{(r)\top}{\bm{K}}^{(r)}/T$ under the normalization
${\bm{K}}^{(r)\top}{\bm{K}}^{(r)}/T=\mathbb{I}_{R}$), the optimization problem
is converted to minimizing the objective function:

$\displaystyle\operatorname{Tr}{\left({\bm{M}}^{(r)\top}{\bm{M}}^{(r)}\right)}-T^{-1}\operatorname{Tr}{\left({\bm{K}}^{(r)\top}{\bm{M}}^{(r)}{\bm{M}}^{(r)\top}{\bm{K}}^{(r)}\right)}.$

Thus, the original local weighted least squares problem is equivalent to
maximizing

$\displaystyle\operatorname{Tr}{\left({\bm{K}}^{(r)\top}{\bm{M}}^{(r)}{\bm{M}}^{(r)\top}{\bm{K}}^{(r)}\right)},$

subject to the restriction
${\bm{K}}^{(r)\top}{\bm{K}}^{(r)}/T=\mathbb{I}_{R}$, which is equivalent to
the optimization problem of the classical factor model.

Our objective is to obtain estimators of the factors and factor loadings. A
two-stage estimation procedure is used to estimate those parameters. Let
$\widehat{{\bm{K}}}^{(r)}$ denote the estimated factor matrix of
${\bm{K}}^{(r)}$, and
$\widehat{{\bm{B}}}_{r}=\left(\widehat{{\bm{b}}}_{1,r},\widehat{{\bm{b}}}_{2,r},\dots,\widehat{{\bm{b}}}_{N,r}\right)^{\top}$
denote the estimator of the time-varying factor loading matrix ${\bm{B}}_{r}$.
Then, $\widehat{{\bm{K}}}^{(r)}$ is $\sqrt{T}$ times eigenvectors
corresponding to the largest $R$ eigenvalues of the $T\times T$ matrix
${\bm{M}}^{(r)}{\bm{M}}^{(r)\top}$, and $\widehat{{\bm{B}}}_{r}$ is
${\bm{M}}^{(r)\top}\widehat{{\bm{K}}}^{(r)}\left(\widehat{{\bm{K}}}^{(r)\top}\widehat{{\bm{K}}}^{(r)}\right)^{-1}$
(it is ${\bm{M}}^{(r)\top}\widehat{{\bm{K}}}^{(r)}/T$ under the condition
${\bm{K}}^{(r)\top}{\bm{K}}^{(r)}/T=\mathbb{I}_{R}$). Therefore, in the first
step, we can acquire estimators $\widehat{{\bm{B}}}_{r}$ of the factor
loadings for $r=1,\dots,T$.

Based on the estimator $\widehat{{\bm{B}}}_{r}$ of the factor loading matrix
obtained in the first stage, we consider another least squares problem in the
second stage to obtain the estimator of the factor ${\bm{k}}_{t}$. The
objective function we would like to minimize is as follows:

$\displaystyle\sum_{x=1}^{N}\left(\ln(m_{x,t})-a_{x}-\widehat{{\bm{b}}}_{x,t}^{\top}{\bm{k}}_{t}\right)^{2}\
\text{for}\ t=1,\dots,T.$

Since we already have $\widehat{{\bm{b}}}_{x,t}$ in the first stage, the
answer to this minimization problem is

$\displaystyle\widehat{{\bm{k}}}_{t}=\left(\sum_{x=1}^{N}\widehat{{\bm{b}}}_{x,t}\widehat{{\bm{b}}}_{x,t}^{\top}\right)^{-1}\left(\sum_{x=1}^{N}\widehat{{\bm{b}}}_{x,t}\left(\ln(m_{x,t})-a_{x}\right)\right)\
\text{for}\ t=1,\dots,T.$

Thus, using the two-stage estimation method, we can obtain consistent
estimators for both the factors and time-varying factor loadings.

Next, we discuss some issues in the kernel estimation.

######  Remark 1.

Boundary kernel. Usually, there exists a boundary bias issue in the kernel
estimation. Instead of using the ordinary kernel function, it is suggested
that a boundary kernel should be used to help us obtain some uniform results.
Let $\lfloor a\rfloor$ represent the greatest integer less than or equal to
$a$, then the boundary kernel we choose to use is as follows:

$\displaystyle
K_{h,tr}^{*}=h^{-1}K_{r}^{*}\left(\frac{t-r}{Th}\right)=\begin{cases}\frac{h^{-1}K\left(\frac{t-r}{Th}\right)}{\int_{-\frac{r}{Th}}^{1}K(u)\mathrm{d}u}&\quad
r\in[1,\lfloor Th\rfloor]\\\ h^{-1}K\left(\frac{t-r}{Th}\right),&\quad
r\in(\lfloor Th\rfloor,T-\lfloor Th\rfloor]\\\
\frac{h^{-1}K\left(\frac{t-r}{Th}\right)}{\int_{-1}^{(1-r/T)/h}K(u)\mathrm{d}u}&\quad
r\in(T-\lfloor Th\rfloor,T]\end{cases}.$

######  Remark 2.

The choice of bandwidth. For the nonparametric local smoothing method, it is
important to determine the bandwidth for the kernel estimation. There are two
ways to choose the bandwidth. The first one is to use a data-driven method,
such as the cross-validation. The other one is to use Silverman’s rule of
thumb to set the bandwidth, which is much easier to compute. Su and Wang
(2017) have shown that choices of the kernel function and the bandwidth have
little impact on the performance of the information criteria. Thus, in the
following empirical analysis, we decide to use the Epanechnikov kernel and its
corresponding Silverman’s rule of thumb bandwidth, which is
$h=(2.35/\sqrt{12})T^{-1/5}N^{-1/10}$.

######  Remark 3.

Determination of the number of factors. There are mainly two methods to
determine the number of factors, $R$. The first one is to use a BIC-type
information criterion proposed by Su and Wang (2017). Under certain
assumptions, the new information criterion can correctly choose the true value
of $R$. However, those assumptions may not hold in real data. Additionally, it
is not easy to implement the out-of-sample forecasting if the chosen value of
$R$ is too large.

The second method is based on the fact that the original local weighted least
squares problem can be transformed into an optimization problem of the
classical factor model. Therefore, the cumulative sum of eigenvalues can help
us identify the number of factors. Let $c$ denote a cut-off value between $0$
and $1$, and $\lambda_{k}$ denote the $k^{th}$ largest eigenvalue of the
matrix ${\bm{M}}^{(r)}{\bm{M}}^{(r)\top}$, then we can choose the value of $R$
as
$\min\\{R:\left(\sum_{k=1}^{R}\lambda_{k}\right)/\left(\sum_{k=1}^{N}\lambda_{k}\right)\geq
c\\}$. In the following analysis, we will set the cut-off value as $c=0.9$ and
empirical analysis shows that only one factor is enough to capture most
characteristics of the mortality data, which is consistent with the Lee and
Carter (1992) model.

###  2.3 Forecasting Method

We now consider how to make out-of-sample forecasting using the time-varying
factor model. Since the factor loadings change over time, we should not only
make predictions of the common factors, but also extrapolate the factor
loadings for each age. We describe the forecasting method for a single factor
model ($R=1$) for the simplicity of notations in the following analysis.
Assume that based on the historical data we have acquired the estimated common
factor and factor loadings using the method mentioned in Section 2.2.

In order to forecast the common factor, we first fit the common factor with an
ARIMA model. Since Akaike Information Criterion (AIC) is asymptotically
equivalent to the cross-validation when the maximum likelihood estimation is
used to fit the model (Stone, 1977), we choose AIC as the model selection
criteria to find the most appropriate ARIMA model. After that, we can use the
chosen model to forecast and obtain prediction intervals for the latent factor
(see more details in Chapter 5 & 9 of Brockwell et al. (1991)).

The factor loading $b_{x,t}$ is assumed to be an unknown piece-wise smooth
function of time $t$. For the purpose of extrapolating the factor loading
$b_{x,t}$ into the future, we will adopt two different methods to achieve the
goal:

  1. 1.

The naive method. We simply assume that in the forecasting horizon,
$b_{x,t}(t>T)$ is set as $b_{x,T}$. This is essentially a parametric
forecasting method and similar to that in the Lee-Carter model but with a
different estimated value. The naive method using constant factor loading has
a simple structure and could provide more stable forecasts in the long term.

  2. 2.

The local regression method. This method is based on a nonparametric
regression method – the local linear regression, to flexibly estimate the
deterministic function $b_{x,t}$ (See more details in Fan and Gijbels (1996)
and Friedman et al. (2001)). Similar method has also been applied in Li and
O’Hare (2017); Li et al. (2015). The local linear regression can easily extend
the most recent trends, which is more suitable for short-term forecasting.

We briefly describe the local regression method in the rest of this section.
The main idea of the local linear regression is to fit the linear regression
using only the observations in the neighbourhood of a target point $t_{0}$.
This so-called localization is achieved by using a weight function
$K_{\lambda}(t,t_{0})=K\left((t-t_{0})/\lambda\right)$, where $K$ is a kernel
function and the index $\lambda$ indicates the width of the neighborhood. One
of the commonly used kernel functions with compact support is Epanechnikov
kernel, which is adopted in this paper. For the Epanechnikov kernel, the
window width parameter $\lambda$ is the radius of the support region, which
can be estimated using out-of-sample validation. The weight function assigns a
weight to each time point $t$ based on the corresponding distance from $t_{0}$
(i.e., $|t-t_{0}|$). In this way, the resulting estimated function is a smooth
function.

Specifically for the forecasting of the time-varying factor loading of each
age $x$, the local linear regression solves a separate weighted least square
problem at each target point $T+h$ ($h=1,2,\dots)$:

$\displaystyle\min_{\alpha(T+h),\beta(T+h)}\sum_{t=1}^{T+h-1}K_{\lambda}(t,T+h)\left(b_{x,t}-\alpha(T+h)-\beta(T+h)t\right)^{2}.$

Note that the notations $\alpha(T+h)$ and $\beta(T+h)$ indicate that the two
parameters under study vary with the point $T+h$ in the local linear method.

Let ${\bm{b}}_{x}=\left(b_{x,1},b_{x,2},\dots,b_{x,T+h-1}\right)^{\top}$,
${\bm{X}}=\begin{pmatrix}1&1&\dots&1\\\ 1&2&\dots&T+h-1\end{pmatrix}^{\top}$,
and ${\bm{W}}(T+h)$ denote the $(T+h-1)\times(T+h-1)$ diagonal matrix with the
$t^{th}$ diagonal element $K_{\lambda}(t,T+h)$. Then by using the weighted
least squares estimation, we can obtain the estimators for $\alpha(T+h)$ and
$\beta(T+h)$ as follows:

$\displaystyle\left(\widehat{\alpha}(T+h),\widehat{\beta}(T+h)\right)^{\top}=\left({\bm{X}}^{\top}{\bm{W}}(T+h){\bm{X}}\right)^{-1}{\bm{X}}^{\top}{\bm{W}}(T+h){\bm{b}}_{x}.$

To ensure that ${\bm{X}}^{\top}{\bm{W}}(T+h){\bm{X}}$ is nonsingular, the
bandwidth parameter $\lambda$ in the kernel function should be selected
properly in practice, see more details in Fan and Gijbels (1996). Therefore
the forecasted factor loading at point $T+h$ is

$\displaystyle\widehat{b}_{x,T+h}=\begin{pmatrix}1&T+h\end{pmatrix}\begin{pmatrix}\widehat{\alpha}(T+h)\\\
\widehat{\beta}(T+h)\end{pmatrix}.$

Note that for $h>1$, the forecasts
$\widehat{b}_{x,T+1},\dots,\widehat{b}_{x,T+h-1}$ are evolved in
${\bm{b}}_{x}$ when estimating the factor loading at the time $T+h$. Following
this method, we can estimate the factor loadings for each age as a smooth
function of time $t$ and then extrapolate the factor loadings into the future.
Combining with the predicted common factors, we can make out-of-sample
predictions of the central death rates using the time-varying factor model.

##  3 Optimal ‘Boundary’ Estimation

Under the framework of time-varying factor model, we assume the factor loading
$b_{x,t}$ is a function of time $t$. In Section 2.3, we introduced two
different methods to extrapolate the factor loading. One is a naive method,
which is more suitable for long-term forecasting; and the other is based on
local linear regression, which is more suitable for short-term forecasting.
Then can we estimate the ‘boundary’ between short-term and long-term
forecasting that divides the forecasting horizon according to the predictive
power of the local regression method and naive method?

To solve this problem, we first propose a new forecasting method, which is a
hybrid of two previously introduced methods. Assume $T_{0}$ is the number of
years used in fitting the model and $k$ ($k=0,1,2,\cdots$) is the optimal
boundary between short-term and long-term forecasting, favoured by the time-
varying models based on the local regression and naive method, respectively.
We have the point forecast estimation of mortality rate $ln(m_{x,t})$ for any
given $x$, $t$ and $k\ (k\geq 1)$ using the hybrid method as

$ln(\widehat{m}_{x,t})=\begin{cases}\widehat{a}_{x}+\widehat{b}_{x,t}\cdot\widehat{k}_{t}&T_{0}+1\leq
t\leq T_{0}+k,\\\
\widehat{a}_{x}+\widehat{b}_{x,T_{0}+k}\cdot\widehat{k}_{t}&t\geq
T_{0}+k+1.\end{cases}$

If $T_{0}+1\leq t\leq T_{0}+k$, the forecast of $\ln(m_{x,t})$ at time $t$ is
$\widehat{a}_{x}+\widehat{b}_{x,t}\cdot\widehat{k}_{t}$, where
$\widehat{b}_{x,t}$ is the extrapolated factor loading at time $t$ based on
the local regression method. When $t\geq T_{0}+k+1$, the forecast at time $t$
is $\widehat{a}_{x}+\widehat{b}_{x,T_{0}+k}\cdot\widehat{k}_{t}$, where
$\widehat{b}_{x,T_{0}+k}$ is time-invariant and obtained using the
extrapolated factor loading at time $T_{0}+k$ based on the local regression
method. For $k=0$, the forecast at time $t\ (t>T_{0})$ is just
$\widehat{a}_{x}+\widehat{b}_{x,T_{0}}\cdot\widehat{k}_{t}$ using the
estimated factor loading at time $T_{0}$. In this case, the hybrid method
degenerates to the naive method. In view of this, $T_{0}+k$ is the time
boundary between short-term and long-term forecasting, and between choosing
the local regression and naive method. Given the value of $k$, the hybrid
method applies the local linear regression for the first $k$ periods in the
forecasting horizon and keeps the factor loadings ($\widehat{b}_{x,T_{0}+k}$)
unchanged thereafter, which combines the local regression and naive methods.
Additionally, the hybrid method guarantees a consistent and smooth transition
from short-term to long-term forecasting.

As discussed in Section 1, different forecasting horizons may favour different
models. Generally, long-term forecasting benefits more from the historical
long-term trend and short-term forecasting relies on the recent trend (Booth
et al., 2002). Since the local linear regression can easily extend the most
recent trend, it is more suitable for short-term forecasting. However, as time
goes by, the recent trend becomes less and less reliable, which is not
suitable for long-term forecasting. On the other hand, the naive method using
constant factor loading is more suitable for long-term forecasting, as it has
a simple structure and would provide more stable forecasts in the long term.
Compared to the classical factor model, the naive method provides more
accurate estimations not only for the factor loadings but also for the common
factors, which helps generate more accurate long-term forecasts.

Based on the hybrid forecasting method, we propose an estimation method of the
optimal ‘boundary’ inspired by Bai (2010). Assume the entire dataset has $T$
years and we consider the first $T_{0}$ years of data as the training set, and
the remaining data with size $T-T_{0}$ as the validation set. Given the value
of $k$, we first fit the time-varying factor model using the training set, and
then apply the hybrid forecasting method to the validation set. We consider
all possible lengths of short-term (long-term) forecasting horizons (i.e. $k$)
and find out an optimal one using least squares estimation. We describe the
estimation procedure as follows.

For the given $x$ and $k$ such that $1\leq k\leq T-T_{0}-1$, define
$\widehat{y}_{x,t}(k)=\widehat{a}_{x}+\widehat{b}_{x,t}\cdot\widehat{k}_{t}$
as the predicted value of $\ln(m_{x,t})$ using the hybrid forecasting method.
When $T_{0}+1\leq t\leq T_{0}+k$, $\widehat{b}_{x,t}$ is forecasted by the
local regression method; And when $T_{0}+k+1\leq t\leq T$,
$\widehat{b}_{x,t}=\widehat{b}_{x,T_{0}+k}$, where $\widehat{b}_{x,T_{0}+k}$
is the predicted factor loading at time $T_{0}+k$ obtained via the local
regression method. Then we define the sum of squared residuals for age $x$ as

$\displaystyle S_{x,T}(k)$
$\displaystyle=\sum_{t=T_{0}+1}^{T}\left(\ln(m_{x,t})-\widehat{y}_{x,t}(k)\right)^{2}$
$\displaystyle=\sum_{t=T_{0}+1}^{T_{0}+k}\left(\ln(m_{x,t})-\widehat{a}_{x}-\widehat{b}_{x,t}\cdot\widehat{k}_{t}\right)^{2}+\sum_{t=T_{0}+k+1}^{T}\left(\ln(m_{x,t})-\widehat{a}_{x}-\widehat{b}_{x,T_{0}+k}\cdot\widehat{k}_{t}\right)^{2},$

where $k=1,2,\dots,T-T_{0}-1$. Here $k$ represents the length of the short-
term forecasting horizon or the ‘boundary’ between short-term (based on the
local regression method) and long-term (based on the naive method)
forecasting. The local linear regression is used to make forecasts from
$T_{0}+1$ to $T_{0}+k$; while the naive method (i.e. assuming
$\widehat{b}_{x,t}$ does not change over the remaining period) is used to make
forecasts from $T_{0}+k+1$ to $T$. We define

$\displaystyle
S_{x,T}(0)=\sum_{t=T_{0}+1}^{T}\left(\ln(m_{x,t})-\widehat{a}_{x}-\widehat{b}_{x,T_{0}}\cdot\widehat{k}_{t}\right)^{2},\
\ \text{for}\ \ k=0$

and

$\displaystyle
S_{x,T}(T-T_{0})=\sum_{t=T_{0}+1}^{T}\left(\ln(m_{x,t})-\widehat{a}_{x}-\widehat{b}_{x,t}\cdot\widehat{k}_{t}\right)^{2},\
\ \text{for}\ \ k=T-T_{0}.$

In this way, $S_{x,T}(k)$ is defined for each $k=0,1,\dots,T-T_{0}$. Thus, the
total sum of squared residuals (SSR) across all ages is defined as

$\displaystyle SSR(k)=\sum_{x=1}^{N}S_{x,T}(k),$

Hence the least squares estimator of the optimal ‘boundary’ is

$\displaystyle\widehat{k}=\underset{0\leq k\leq T-T_{0}}{argmin}SSR(k).$

The estimated optimal ‘boundary’ between the short-term (based on local linear
regression) and long-term (based on naive method) forecasting is the time
$\widehat{k}$ that leads to the smallest SSR.

##  4 Data

The mortality data used in this paper are extracted from the Human Mortality
Database (HMD) (37). Six countries are selected for the empirical analysis in
Section 4 and Section 5. For each country, age-sex-specific death rates are
available annually for the entire population. The selected countries are shown
in Table 1 along with the corresponding available time horizons, which will be
used for empirical analysis.

Table 1: Time horizon for different countries Country |  start year |  end year |  length
---|---|---|---
AUSTRALIA |  1921 |  2018 |  98
CANADA |  1921 |  2016 |  96
FRANCE |  1816 |  2017 |  202
ITALY |  1872 |  2017 |  146
JAPAN |  1947 |  2018 |  72
USA |  1933 |  2017 |  85

The mortality data are generally available from age $0$ to age $110+$ for each
year. Since measures of mortality at very old ages are unreliable (Lee and
Carter, 1992), we decide not to use mortality data of age $91$ and over in the
following analysis and end up with $N=91$ ages.

In order to investigate whether the factor loadings are time-varying or time-
invariant in the empirical data. We conduct an exploratory data analysis by
applying the Lee-Carter Model on the US mortality data with rolling-window
time frames. We first divide the entire dataset into 44 subsets (each with 40
yearly observations) with the first subset from year 1933 to year 1972, the
second subset from year 1934 to year 1973, and so on. We then fit the Lee-
Carter model on each of the subset and extract the factor loading $b_{x}$ for
each time frame. We plot the factor loadings of some selected ages in Figure
1. We can see that the factor loadings possess different dynamic patterns for
different ages and they are not time-invariant.

Figure 1: Factor loadings for ages 20, 40, 60, 80 over 44 rolling-window time
frames

##  5 Empirical Results and Analysis

In the first two subsections, we present the application results of the time-
varying factor model using age-specific mortality data of the US. We compare
the time-varying factor models based on both the naive and local regression
forecasting methods with Lee-Carter model (the classical factor model with one
factor) via out-of-sample forecasting performance. Empirical results by gender
are provided in Appendix A. Section 5.3 further compares the forecasting
performance across multiple countries and models based on different fitting
and forecasting horizons. And in the Section 5.4, we estimate the optimal
‘boundary’ between short-term and long-term forecasting for different
countries.

###  5.1 Model Fitting

We fit the US mortality data using the estimation method of the time-varying
factor model introduced in Section 2. The number of factors estimated is 1
($\widehat{R}=1$), which is consistent with the Lee-Carter model. More
empirical results of the time-varying model with multiple factors are shown in
Appendix C.

Using the model selection criteria AIC, we find that the common factor
$k_{t}$, obtained from the time-varying factor model, follows an ARIMA
$(1,1,0)$ with drift model. This model can capture most of the characteristics
of the common factor. Our fitted model of the common factor $k_{t}$ is as
follows:

$\displaystyle\bigtriangledown
k_{t}=\underset{(0.2791)}{-1.4116}+\underset{(0.1046)}{0.3271}\bigtriangledown
k_{t-1}+e_{t},$

where $\bigtriangledown$ refers to the first order differencing and $e_{t}$
represents the error term. The numbers in the parentheses are the standard
errors of the corresponding parameters. With the ARIMA model built above, we
can then forecast the common factor into the future.

As a comparison, we also list the fitted ARIMA model of the common factor (The
number of factors estimated is 1.) using the classical factor model111In this
case, the classical factor model has the same model structure as the Lee-
Carter model with $R=1$. below:

$\displaystyle\bigtriangledown
k_{t}=\underset{(0.2837)}{-1.4046}+\underset{(0.1051)}{0.3114}\bigtriangledown
k_{t-1}+e_{t}.$

We can see that the ARIMA models of the common factors estimated from the
time-varying factor model and the Lee-Carter model are close. The estimated
common factors, plotted in Figure 2, tend to decrease linearly and show
similar dynamic patterns. The common factor is regarded as the index of the
level of mortality, which captures major influence on death rates of all ages.

Figure 2: Plots of the estimated common factors for the time-varying factor
model & the Lee-Carter model

Figure 4 displays the comparison of the factor loadings between the time-
varying factor model and the Lee-Carter model for selected ages. Compared with
time-invariant factor loadings (the dashed lines), the time-varying factor
loadings (the solid curves) change smoothly overtime, see Figure 4. It is
interesting to notice that, no matter which age it is, the corresponding
factor loadings always reach their own minimum or maximum values during
$1960$s or $1970$s. For older people (over age $40$), the factor loadings
usually arrive at their maximum values during $1960$s or $1970$s, which means
the death rates of older people are more sensitive to the latent common factor
during that period. For the younger ages (below age $40$), however, the
corresponding factor loadings reach their minimum values during the same
period, which means the death rates of younger people are less sensitive to
the latent factor during that time. The only exception is the factor loadings
of the infant group, whose dynamic pattern is more similar to that of the
older group.

Figure 3: Plots of the estimated time-invariant factor loadings (dashed lines)
& the time-varying factor loadings (solid lines) for age $0,10,\dots,90$.

Figure 4: The actual data (black solid lines) versus the fitted values from
the time-varying model (red dashed lines) and the Lee-Carter model (blue
dotted lines); the data have been log-transformed & demeaned.

Figure 4 shows the fitted death rates of both the time-varying factor model
and the Lee-Carter model with empirical observations for selected ages.
Obviously, no matter which age it is, the time-varying factor model fits
better than the Lee-Carter model. We use the mean squared error (MSE)222The
MSE for the time-varying model is computed as follows:
$\displaystyle\text{MSE}=\frac{1}{NT}\sum_{x}\sum_{t}\left(\ln(m_{x,t})-a_{x}-\widehat{{\bm{b}}}_{x,t}^{\top}\widehat{{\bm{k}}}_{t}\right)^{2}$
Computation of the MSE for the Lee-Carter model is the same as above except
that $\widehat{{\bm{b}}}_{x,t}$ is replaced by $\widehat{{\bm{b}}}_{x}$. to
evaluate the goodness of fit. As a result, the overall MSE of the time-varying
factor model is $0.001990$, which is much smaller than that of the Lee-Carter
model, $0.006690$ (three times bigger than the former one). Therefore, the
time-varying factor model performs much better than the Lee-Carter with
respect to the in-sample fitting.

Although the time-varying factor model works better in the fitting procedure,
the problem of overfitting may exist due to the increased complexity of the
model. Through the Monte Carlo simulation studies in Section 6, we will see
that overfitting is harmful to forecasting. Usually, an overfitting model
performs too well in the fitting sample to have good generalization ability in
forecasting. Generally speaking, we can always improve a model’s in-sample
fitting performance by increasing the complexity of the model, which, however,
cannot guarantee a better forecasting performance in the future. Thus we will
use the out-of-sample validation method to investigate whether the time-
varying model can enhance the out-of-sample forecasting performance in the
next subsection.

###  5.2 Out-of-sample Forecasting

In this subsection, we use the original US mortality data over the first 60
years as the training set (from 1933 to 1992) to fit the models, and then
forecast the mortality rates in the testing set (from 1993 to 2017) using the
fitted models. The predicted values are compared with the actual data in the
testing set to see which model is better at the out-of-sample forecasting. We
apply the mean squared prediction error (MSPE) as the measure to evaluate the
out-of-sample forecasting performance.

Figure 5: Out-of-sample forecast of the common factor, with the model fitted
on $1933$ to $1999$ and the forecast horizon over $1993$ to $2017$; predicted
value (dashed line), $80\%$ PI (dotted line)

Figure 5 plots the historical and predicted values of the common factor of the
time varying factor model along with the associated $80\%$ prediction
intervals, which is based on the ARIMA model fitted in Section 5.1. The dashed
downward line shows that the latent factor will keep declining in the future,
and there is an $80\%$ chance that a future observation will be covered by the
corresponding prediction interval (represented by area between the red dashed
lines).

Since the time-varying factor model and the Lee-Carter model have similar
common factor and the corresponding fitted ARIMA model, the forecasts of the
common factor are close to each other too. Hence, the major difference of
prediction accuracy between the time-varying factor model and the Lee-Carter
model lies in the factor loadings. We extrapolate the factor loadings obtained
from the time-varying factor model using both the naive method and local
linear regression introduced in Section 2.3, respectively. We then forecast
the mortality rates into the future using both the time-varying and Lee-Carter
models.

Figure 6: Plots of the estimated and extrapolated factor loadings based on
naive method (dotted lines) & local regression method (dashed lines) for age
$0,10,\dots,90$.

Figure 7: The actual data (solid line) versus the predicted values from the
time-varying model (naive method: dash-dotted line; local linear regression:
dotted line) and the Lee-Carter model (dashed line); the data have been log-
transformed.

Figure 7 plots the estimated and extrapolated factor loadings of the time
varying factor models. Figure 7 plots the actual data and predicted values
using the three above-mentioned methods. From Figure 7, we see that the local
regression method follows the recent historical trend of factor loadings,
while the naive method stays at a constant level. Theoretically, if future
factor loadings do not deviate significantly from the recent historical trend,
the local linear regression may perform better than the naive method in
forecasting. However, it may only be reasonable to assume that factor loadings
will follow the local trends in the short-term. For long-term forecasting,
this assumption is less reliable and the non-parametric forecasting method
would lead to inferior results. In Section 5.3, we observe similar results in
other countries’ mortality forecasting. Hence, the naive method, with time-
invariant forecasted factor loadings, is more suitable for long-term
forecasting. And it may also be suitable for short-term forecasting if the
long-term trend is consistent with the short-term trend. Since the benefit of
using the local regression method decreases as the forecasting horizon
increases, it is worthwhile to ask whether an optimal forecasting horizon
exists for using the local regression method. This question will be answered
in Section 5.4.

The empirical analysis suggests that, the time-varying factor model (based on
the naive method) performs better than the Lee-Carter model over the entire
forecasting horizon ($1993$-$2017$). Using the mean squared prediction error
(MSPE) to evaluate the out-of-sample forecasting performance, we see that the
overall MSPE for the Lee-Carter model is $0.03085$, while the overall MSPE for
the time-varying factor model (based on the naive method) is only $0.01804$.
However, if we choose to use the local linear regression to extrapolate factor
loadings, the time-varying factor model performs worse than Lee-Carter model,
with the MSPE being $0.04768$.

Figure 8 shows the year-specific MSPE for the time-varying factor models and
the Lee-Carter model over the forecasting horizon $1993$ to $2017$. The year-
specific MSPE is computed by averaging MSPE over all ages for each forecasting
year. From Figure 8, we can see that for the majority years, the MSPE of the
time-varying model with the naive forecasting method is always the smallest
one. From Section 5.4, we can see the reason is that the optimal ’boundary’
between short-term and long-term forecasting in this case is estimated to be
$0$. So the the time-varying model with the naive forecasting method is the
best for both short-term and long-term forecasting. The MSPE for all the three
methods are generally increasing over the years, as it is harder to forecast
the farther future. We also notice that the time-varying factor model based on
the local regression method works better than the Lee-Carter model over 1993
to 1995. However, it has the worst performance for long-term forecasting. The
time-varying model based on local regression method assumes the factor
loadings change over time in the future, but it can only extend the recent
trend, which may not be suitable for long-term forecasting. On the other hand,
the Lee-Carter model and the time-varying model based on the naive method
extrapolate factor loadings into the future as constants, which are usually
more suitable for long-term forecasting.

Figure 8: US: Year-specific MSPE for the time-varying model and the Lee-Carter
model over $1993$ to $2017$; for time-varying model, both the naive method and
the local regression method are used Figure 9: US: Age-specific MSPE for both
the time-varying model and the Lee-Carter model; for time-varying model, both
naive method and local regression method are used

Next, we investigate the forecasting performance of the time-varying factor
model at different ages. Figure 9 shows the age-specific MSPE for the time-
varying factor models and the Lee-Carter model. The square-dashed lines and
triangle-dotted lines represent the MSPE of the time-varying factor models
with the naive and local regression methods for each age, respectively, and
the circle-solid line represents the MSPE of the Lee-Carter model. We find
that the naive forecasting method based on the time-varying factor model is
almost always better than the local regression method for any age in this
data. And roughly speaking, no matter which extrapolation method we choose to
use, the time-varying models provide more accurate forecasts than the Lee-
Carter model for age groups $20\sim 40$ and $60\sim 80$. However, for the age
group $40\sim 60$, the forecasting performance of the time-varying factor
models is worse than that of the Lee-Carter model. Thus, even though by using
naive extrapolation method the time-varying factor model improves the overall
performance (over $40\%$ in terms of the MSPE) significantly, it cannot
outperform the Lee-Carter model for some ages. The main advantage of the time-
varying model is to forecast mortality rates for the young adulthood ($20\sim
40$) and the older adulthood ($60\sim 80$).

###  5.3 Model Comparisons for Multiple Countries

We apply and compare different models using mortality data of multiple
countries. In particular, the functional data model proposed by Hyndman and
Ullah (2007) is also considered for comparison purposes and we call it
‘Hyndman-Ullah model’ in the rest of the paper. It is a multi-factor extension
of the Lee-Carter model, allowing for multiple age-time interaction terms to
capture the complex structure of the data. Similar to the Lee-Carter model
(i.e. the classical factor model with one factor), it is also commonly
considered as a benchmark for mortality forecasting. To determine the order
$K$ of Hyndman-Ullah model (similar to the number of factors $R$ in this
paper), we choose the best value of $K$ by minimizing Integrated Squared
Forest Error. Please refer to Hyndman and Ullah (2007) for more details.

In Table LABEL:t2, we present the results of the overall MSPE for different
countries, forecast horizons and methods. We use the longest available dataset
for training purposes, with different forecasting horizons listed in Table
LABEL:t2. To investigate the performance of the short-term and long-term
forecasting, we consider multiple forecasting horizons with different lengths,
including $5$, $10$, $15$, $20$ and $25$ years. Note that, as the number of
total historical years is fixed for each country, the number of training years
changes with the length of the forecasting horizon. For example, the total
historical years for Australia are 1921$\sim$2018 (98 years). Therefore, the
corresponding training years for forecasting horizons 2014$\sim$2018 and
2009$\sim$2018 are 1921$\sim$2013 (93 years) and 1921$\sim$2008 (88 years),
respectively. Please refer to Table 1 for the available historical years for
each country. The number of factors is always estimated as 1 ($\widehat{R}=1$)
for both the classical factor model the time-varying model. Hence, we use
‘Lee-Carter’ to represent the classical factor model (with $R=1$) in this
section. As for Hyndman-Ullah model, we record the estimated order $K$ for
each scenario in Table LABEL:t2.

t2

From Table LABEL:t2, we see that the time-varying models can significantly
improve the out-of-sample forecasting performance compared with the Lee-Carter
model. The time-varying factor models and Hyndman-Ullah model perform the best
in most cases. The Hyndman-Ullah model is especially suitable for mortality
forecasting of the France data, as shown in Hyndman and Ullah (2007). In
addition, for the time-varying factor models, the local regression method
tends to perform better for short-term forecasting, while the naive method is
better for long-term forecasting. Compared with all the other models, the
time-varying model based on the naive method show superior results for long-
term forecasting, while the Lee-Carter model always has the worst forecasting
performance.

Figure 10: Year-specific MSPE by country and method (Hyndman-Ullah: functional
data model in Hyndman and Ullah (2007), Lee-Carter: classical factor model,
local regression: time-varying factor model based on local linear regression,
naive: time-varying factor model based on naive method); length of forecast
horizon is 25 years

In Figure 10, by fixing the length of the forecast horizon to be 25 years, we
plot the year-specific MSPE for all the countries using different forecasting
methods. Roughly speaking, the time-varying model based on the naive method
always has the most accurate forecasts, while the Lee-Carter model always
performs the worst. Additionally, the time-varying model based on the local
regression method usually performs well (similar to or better than the naive
method) at the first few years, while it deteriorates as time goes by. For
mortality data of Australia and the US, it has worse forecasting performances
than the Lee-Carter model in the long term.

###  5.4 Estimate the Optimal ‘Boundary’

From the empirical results above, we see that under the framework of time-
varying model, the local regression method (by assuming factor loading will
change in the future) is better at short-term forecasting while the naive
method (by assuming factor loading is in-variant in the future) is better at
long-term forecasting. Therefore, we are interested in the boundary between
short-term and long-term forecasting that divides the forecast horizon
according to the predictive power of the local regression method and the naive
method.

By applying the estimation method introduced in Section 3, we estimate the
optimal ‘boundary’ between the short-term and long-term forecasting, which is
favoured by the local regression method (time-varying forecast of the factor
loading) and the naive method (time-invariant forecast of the factor loading)
respectively. Recall that the optimal ‘boundary’ can be regarded as the
optimal value of $k$ defined in the hybrid forecasting method in Section 3.
Applying the same datasets introduced in Section 4, we use the last 25 years
of the historical data as the validation set, and the remaining data as the
training set. In addition, to check the sensitivity of the least squares
estimator to the division of the validation and testing sets, we consider the
last $p$ years ($p=15,20,25,30$, respectively) of data as the validation set
and the remaining data as the training set. The results are shown in Appendix
B.

Figure 11: Plots of the total sum of squared residuals (SSR) versus the length
($k$) of the short-term forecast horizon (based on the hybrid forecasting
method of time-varying factor model); length of forecast horizon: $25$

As defined in Section 3, for each set of data we compute the least squares
estimator for the optimal ‘boundary’ in the forecasting horizon as

$\displaystyle\widehat{k}=\underset{0\leq k\leq T-T_{0}}{argmin}SSR(k),$

where $k=0,1,2,\dots,T-T_{0}$. Here $k$ represents the length of the short-
term forecasting horizon, or the ‘boundary’ between short-term (based on the
local regression method) and long-term (based on the naive method)
forecasting. The local linear regression is used to make forecasts from
$T_{0}+1$ to $T_{0}+k$; while the naive method (i.e. keep the factor loading
the same as that in time $T_{0}+k$) is used to make forecasts from $T_{0}+k+1$
to $T$.

We plot the total SSR versus $k$ in Figure 11. As shown in the plots, for each
country, there is a minimum point of $k$ corresponding to the smallest values
of SSR, which indicates the optimal length of the short-term forecasting
horizon for the time-varying factor model. For example, the plot of Italy
shows that, when the value of $k$ equals 7, SSR achieves the smallest value.
Thus, based on the historical Italy mortality data, we suggest that it is
better to use the local regression method for the short-term forecasting (less
than 7 leads), as it puts more weights on the recent observations. As for
forecasting more than $7$ years ahead, the naive method generates more
accurate predictions. However, as for the US data, the optimal length of the
short-term forecasting horizon is 0, which means there is no need to
extrapolate factor loadings using local linear regression and we should use
the naive method alone. In Figure 11, we can observe jumps at $k=1$ for some
countries (such as Canada and France). When $k=0$, the hybrid method is just
the naive method and the extrapolation of factor loading is based on the
historical estimation; while when $k>1$, the hybrid method is based on the
local regression method before the time $T_{0}+k$. Therefore, a sudden change
of SSR when $k$ increases from $0$ to $1$ indicates that the pattern of factor
loadings has a relatively large change after the time $T_{0}$.

Now we have observed the existence of positive optimal ‘boundary’ for some
countries (such as Australia, Italy and Japan). Then based on the estimation
of optimal ‘boundary’ $k$, we can compare the out-of-sample forecasting
performance of the hybrid method with other methods. Here, we consider the
Australian mortality data as an example since the estimation of the optimal
boundary for Australian data is relatively stable (see details in Remark 4 and
Appendix B).

Similar to previous analysis, we choose the last 25 years of the historical
data as the testing set. And the remaining data is the training set. To
estimate $k$, we use the last 25 years of the training data as the validation
set to estimate the optimal ‘boundary’. We then apply the estimation procedure
in Section 3 to the training data and the forecasting method in the testing
data using the estimated optimal ‘boundary’ to get the out-of-sample
performance. The empirical results are shown in Figure 12 and Table 2.

In Table 2, we compute the overall MSPE over $1994$ to $2018$ for each
forecasting method. We find that the time-varying model based on the hybrid
forecasting method has the best out-of-sample performance among all four
methods. And the naive method is a little bit worse than the hybrid method.
Additionally, in Figure 12, we plot the year-specific MSPE for all different
methods. For short-term forecasting, the hybrid method shows similar
performance with the local regression method; while for the long-term
forecasting, the hybrid method shows similar performance with the naive
method. Therefore, the hybrid method is superior for both short-term and long-
term forecasting, as it benefits from the advantages of both methods. In
practice, both naive method and hybrid method are recommended. The hybrid
method could produce more accurate predictions while the naive method is much
easier to implement.

Table 2: Australia: Overall MSPE over $1994$ to $2018$ (Lee-Carter: classical factor model, TV-Local Regression: time-varying factor model based on local regression method), TV-Naive: time-varying factor model based on naive method), TV-Hybrid: time-varying factor model based on hybrid method) Country |  Lee-Carter |  TV-Local Regression |  TV-Naive |  TV-Hybrid
---|---|---|---|---
Australia |  0.11162 |  0.08932 |  0.04342 |  0.04288
Figure 12: Australia: Year-specific MSPE over $1994$ to $2018$ (Lee-Carter:
classical factor model, local regression: time-varying factor model based on
local regression method), naive: time-varying factor model based on naive
method), hybrid: time-varying factor model based on hybrid method)

######  Remark 4.

In Appendix B, we further investigate the sensitivity of the optimal
‘boundary’ estimation by showing plots of the total sum of squared residuals
(SSR) versus $k$ using different lengths of training and testing datasets. For
some countries, like Japan and USA, the least squares estimators are
relatively stable. However, for countries like Canada and Italy, the
estimation of the optimal boundary is somewhat sensitive to the length of the
dataset. This is reasonable as the time-varying factor model highly depends on
the data and its intrinsic patterns. Therefore, when using different training
data and forecasting horizons to estimate optimal ‘boundary’, we may obtain
different values of $k$.

##  6 Monte Carlo Simulations

In this section, we further investigate the prediction performance of the
time-varying factor model and the classical factor model through Monte Carlo
simulations. We use examples with different structures of the factor loadings
to illustrate that the time-varying factor model can improve the forecasting
accuracy when the ‘true’ factor loadings change over time. In addition, we
explain under which conditions the naive method performs better than the local
regression method even in the short-term.

Similar to the previous empirical analysis, we denote the classical factor
model as ‘Lee-Carter’ since there is only one factor. We denote the two
forecasting methods based on the time-varying factor model as ‘TV-Local
Regression’ method and ‘TV-Naive’ method, respectively. We show that when the
‘true’ factor loadings change over time, both forecasting methods based on the
time-varying factor model outperform the Lee-Carter method in forecasting.
Moreover, the ‘TV-Naive’ method performs similarly to the ‘TV-Local
Regression’ method for short-term forecasting; while it performs better than
the ‘TV-Local Regression’ for long-term forecasting.

###  6.1 Data Generating Processes (DGP’s)

We generate the centered data $x_{i,t}$ with one common factor:

$\displaystyle x_{i,t}=b_{i,t}\cdot{k_{t}}+\epsilon_{i,t},\ \ i=1,2,\ldots,N,$

where $k_{t}=k_{t-1}+w_{t}$. $w_{t}$ follows independent identically
distributed normal distributions, $N(0,0.8^{2})$. Thus the common factor
$k_{t}$ follows a random walk. We consider the following settings for
different factor loadings $b_{i,t}$ and error terms $\epsilon_{i,t}$. In each
setting below, we apply the normalization condition mentioned in Section 2, so
that $b_{i,t}$ is normalized to sum to unity for each $t$.

  * •

DGP 1 (time-invariant factor loading):

$\displaystyle b_{i,t}=b_{i}\sim{i.i.d}\ {uniform(0,1)},\
\epsilon_{i,t}\sim{i.i.d}\ N(0,0.1^{2}).$

  * •

DGP 2 (single-point structural change):

For $i=1,2,\cdots,N/2$,

$\displaystyle b_{i,t}=\left\\{\begin{array}[]{ll}b_{i}&\quad for\
t=1,2,...,T/2\\\ b_{i}+1&\quad for\ t=T/2+1,...,T\end{array}\right.;$

and for $i=N/2+1,\cdots,N$,

$\displaystyle b_{i,t}=\left\\{\begin{array}[]{ll}b_{i}&\quad for\
t=1,2,...,T/2\\\ b_{i}-1&\quad for\ t=T/2+1,...,T\end{array}\right.;$
$\displaystyle b_{i}\sim{i.i.d}\ {uniform(1.1,1.9)},\
\epsilon_{i,t}\sim{i.i.d}\ N(0,0.03^{2}).$

  * •

DGP 3 (continuous structural change):

$\displaystyle b_{i,t}=\frac{1}{1+e^{(\frac{6i}{N}+2-\frac{12t}{T})}},\
\epsilon_{i,t}\sim{i.i.d}\ N(0,0.1^{2}).$

DGP 1 follows the Lee-Carter model with time-invariant factor loadings, and
DGP 2 and DGP 3 exhibit different structures of the time-varying factor
loadings. DGP 2 describes a single-point structural change in the factor
loadings; while DGP 3 considers a continuous structural change in the factor
loadings. For each $i$, the factor loadings generated in DGP 3 are monotonic
functions and would converge to some constant as time goes by. The factor
loadings in DGP 3 may go up and down and would never diverge to extreme
values, which is similar to the estimated time-varying factor loadings in
Figure 4 using the US mortality data.

###  6.2 Comparison of the Forecasting Performance

To compare the forecasting performance of the time-varying factor model and
the Lee-Carter model (classical time-invariant factor model), we use the out-
of-sample testing approach in the following analysis. For each DGP, we
simulate 100 data sets with the dimension and sample sizes $N=T=100$. For each
data set, we consider the first $k$ years of the data as the training set, and
the remaining $T-k$ years of the data as the testing set
($k=70,75,80,85,90,95$, respectively). The model is first fitted using the
training set, and then forecasted in the testing set. We employ the mean
squared prediction error (MSPE) as the measure to evaluate performance of
different models.

Table 3 reports the comparison results based on various lengths of the
training and testing sets. An example of the estimated and forecasted factor
loadings is shown in Figure 13 to better explain the results. As shown in
Table 3, the two time-varying methods perform worse than the Lee-Carter model
for DGP 1, which assumes the time-invariant factor loadings. The less accurate
prediction results can be attributed to the inaccurate estimation from the
time-varying model, which is supported by the left plot of Figure 13. When the
‘true’ factor loadings are time-invariant, the estimation from the time-
varying model goes up and down randomly due to the over-fitting problem of the
non-parametric estimating method. Therefore, the forecasting based on these
estimation is not satisfied. However, the Lee-Carter method provides a close
estimation and better forecasting in this case.

Table 3: Comparison of forecasting performance of the time-varying factor model and the Lee-Carter model (based on the different lengths of training sets) DGPs | methods |  70 |  75 |  80 |  85 |  90 |  95
---|---|---|---|---|---|---|---
| TV-Local Regression |  $0.2300$ |  $0.1873$ |  $0.1249$ |  $0.0852$ |  $0.0590$ |  $0.0344$
DGP 1 | TV-Naive |  $0.2239$ |  $0.1849$ |  $0.1228$ |  $0.0837$ |  $0.0582$ |  $0.0342$
| Lee-Carter |  $\bf{0.2209}$ |  $\bf{0.1839}$ |  $\bf{0.1222}$ |  $\bf{0.0829}$ |  $\bf{0.0575}$ |  $\bf{0.0336}$
| TV-Local Regression |  $0.2597$ |  $0.2046$ |  $0.1230$ |  $0.0818$ |  $0.0528$ |  $0.0265$
DGP 2 | TV-Naive |  $\bf{0.2291}$ |  $\bf{0.1874}$ |  $\bf{0.1227}$ |  $\bf{0.0817}$ |  $\bf{0.0527}$ |  $\bf{0.0265}$
| Lee-Carter |  $0.2542$ |  $0.2121$ |  $0.1482$ |  $0.0946$ |  $0.0643$ |  $0.0372$
| TV-Local Regression |  $0.1757$ |  $0.1384$ |  $0.0987$ |  $0.0737$ |  $0.0506$ |  $0.0365$
DGP 3 | TV-Naive |  $\bf{0.1704}$ |  $\bf{0.1319}$ |  $\bf{0.0940}$ |  $\bf{0.0724}$ |  $\bf{0.0499}$ |  $\bf{0.0362}$
| Lee-Carter |  $0.2078$ |  $0.1759$ |  $0.1456$ |  $0.1162$ |  $0.0917$ |  $0.0748$

Figure 13: Comparison of the factor loadings: estimation and forecast. From
left to right: DGP 1, DGP 2, DGP 3. $k=70$. Solid line: true factor loadings.
Dot-dashed line: estimation from the classical factor model (‘Lee-Carter’).
Dotted line: estimation from the time-varying factor model. Long dashed line:
‘TV-Local Regression’. Short dashed line: ‘TV-Naive’.

DGP 2 and DGP 3 follow the structures of time-varying factor models with
abrupt and continuous changes in the factor loadings, respectively. From Table
3, we see that both the two time-varying methods perform better than the Lee-
Carter method when the ‘true’ factor loadings are time-varying. In particular,
the naive method performs the best in these two cases, especially for long-
term forecasting. The superior forecasts of the two time-varying methods
result from the more accurate estimation of the time-varying factor loadings,
which can be seen from the middle and right plots in Figure 13. We find that
the Lee-Carter model cannot capture the changed factor loadings as it assumes
the factor loadings are time-invariant, while the time-varying model can
estimate those changing factor loadings accurately and provides a solid
foundation for the forecasting step.

Further analysis of Table 3 suggests that the local regression method and the
naive method preform similarly for the relatively short-term forecasting
($k=85,90,95$), while for long-term forecasting ($k=30,25$), the naive method
performs better. This phenomenon can be explained by the plot of DGP 3 in
Figure 13. From the plot, we see that the forecast of the local regression
method follows the trend of the estimated factor loadings. Therefore, when the
trend remains the same in a short period of time, the forecast of the local
regression method is satisfying. However, as the trend of the ‘true’ factor
loadings changes, the forecast of the local regression method diverges away
from the true values in the long term. This is the drawback of the non-
parametric forecasting method. On the other hand, the constant factor loadings
used in the forecasting procedure of the naive method guarantee that the
factor loadings will not diverge dramatically and result in a better
performance for long-term forecasting. In addition, although the ‘true’ factor
loading in the training set changes, if it is time-invariant in the
forecasting horizon, the naive method is also better than the local regression
method even in the short term (which is supported by the plot of DGP 2 in
Figure 13). The aforementioned analysis could be used to explain why the
estimated optimal ‘boundary’ of the US data is 0 in Section 5.4. From Figure 4
we see that in the first years of the forecasting horizon, the trends of the
factor loadings are either different from that in the training set or remain
flat. However, the local regression method cannot capture the unknown changing
trends, so the naive method can outperform it even in the short term. In this
case, the naive method not only captures the time-varying factor loadings in
the estimation step, but also uses the constant factor loadings in the
forecasting step.

##  7 Conclusion

This paper develops a time-varying factor model for mortality modelling. Two
forecasting methods, the local linear regression and naive method, are used to
extrapolate the factor loadings into the future. To understand the optimal
forecasting horizon of the two forecasting methods, we propose an approach to
estimate the optimal ‘boundary’ between short-term and long-term forecasting,
which is favoured by the local linear regression and the naive method,
respectively. Empirical analysis on mortality forecasting of multiple
countries demonstrates the advantages of the new model.

Like other existing mortality models, the time-varying factor model proposed
in this paper can only extract linear features from the original data. For
future research, we are interested in developing a general non-parametric
factor model, which can capture non-linear features in mortality data. Modern
statistical techniques, such as neural networks and non-parametric estimation
methods, can be used to fulfil this target. In addition, we are interested in
extending the age range considered in this research to include more advanced
ages (Huang et al., 2020) and exploring time-varying factor models for multi-
population mortality forecasting.

## References

  * (1)
  * Anderson (1963) Anderson, T. W. (1963), ‘The use of factor analysis in the statistical analysis of multiple time series’, Psychometrika 28(1), 1–25.
  * Bai (2009) Bai, J. (2009), ‘Panel data models with interactive fixed effects’, Econometrica 77(4), 1229–1279.
  * Bai (2010) Bai, J. (2010), ‘Common breaks in means and variances for panel data’, Journal of Econometrics 157(1), 78–92.
  * Bai and Ng (2002) Bai, J. and Ng, S. (2002), ‘Determining the number of factors in approximate factor models’, Econometrica 70(1), 191–221.
  * Bell (1997) Bell, W. R. (1997), ‘Comparing and assessing time series methods for forecasting age-specific fertility and mortality rates’, Journal of Official Statistics 13(3), 279.
  * Booth et al. (2006) Booth, H., Hyndman, R. J., Tickle, L. and De Jong, P. (2006), ‘Lee-Carter mortality forecasting: a multi-country comparison of variants and extensions’, Demographic Research 15, 289–310.
  * Booth et al. (2002) Booth, H., Maindonald, J. and Smith, L. (2002), ‘Applying Lee-Carter under conditions of variable mortality decline’, Population Studies 56(3), 325–336.
  * Booth and Tickle (2008) Booth, H. and Tickle, L. (2008), ‘Mortality modelling and forecasting: A review of methods’, Annals of Actuarial Science 3(1-2), 3–43.
  * Booth et al. (2004) Booth, H., Tickle, L. et al. (2004), ‘Beyond three score years and ten: prospects for longevity in Australia’, People and Place 12(1), 15.
  * Breitung and Eickmeier (2011) Breitung, J. and Eickmeier, S. (2011), ‘Testing for structural breaks in dynamic factor models’, Journal of Econometrics 163(1), 71 – 84.
  * Brockwell et al. (1991) Brockwell, P. J., Davis, R. A. and Fienberg, S. E. (1991), Time Series: Theory and Methods: Theory and Methods, Springer Science & Business Media.
  * Cairns et al. (2009) Cairns, A. J. G., Blake, D., Dowd, K., Coughlan, G. D., Epstein, D., Ong, A. and Balevich, I. (2009), ‘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.
  * Chang et al. (2018) Chang, J., Guo, B., Yao, Q. et al. (2018), ‘Principal component analysis for second-order stationary vector time series’, The Annals of Statistics 46(5), 2094–2124.
  * Chen et al. (2014) Chen, L., Dolado, J. J. and Gonzalo, J. (2014), ‘Detecting big structural breaks in large factor models’, Journal of Econometrics 180(1), 30 – 48.
  * CMI (2009) CMI (2009), Continuous Mortality Investigation: A prototype mortality projections model: part two – detailed analysis, Working Paper 39, The Institute of Actuaries and the Faculty of Actuaries.
  * Currie et al. (2004) Currie, I. D., Durban, M. and Eilers, P. H. (2004), ‘Smoothing and forecasting mortality rates’, Statistical Modelling 4(4), 279–298.
  * Fan and Gijbels (1996) Fan, J. and Gijbels, I. (1996), Local polynomial modelling and its applications: monographs on statistics and applied probability 66, Vol. 66, CRC Press.
  * Friedman et al. (2001) Friedman, J., Hastie, T. and Tibshirani, R. (2001), The Elements of Statistical Learning, Vol. 1, Springer series in statistics New York, NY, USA:.
  * Hollmann et al. (1999) Hollmann, F. W., Mulder, T. J. and Kallan, J. E. (1999), Methodology & Assumptions for the Population Projections of the United States: 1999 to 2010, US Department of Commerce, Bureau of the Census, Population Division ….
  * Huang et al. (2020) Huang, F., Maller, R. and Ning, X. (2020), ‘Modelling life tables with advanced ages: An extreme value theory approach’, Insurance: Mathematics and Economics 93, 95 – 115.
  * Hyndman and Ullah (2007) Hyndman, R. J. and Ullah, M. S. (2007), ‘Robust forecasting of mortality and fertility rates: a functional data approach’, Computational Statistics & Data Analysis 51(10), 4942–4956.
  * Lam and Yao (2012) Lam, C. and Yao, Q. (2012), ‘Factor modeling for high-dimensional time series: inference for the number of factors’, The Annals of Statistics 40(2), 694–726.
  * Lee and Carter (1992) Lee, R. D. and Carter, L. R. (1992), ‘Modeling and forecasting US mortality’, Journal of the American Statistical Association 87(419), 659–671.
  * Lee and Miller (2001) Lee, R. and Miller, T. (2001), ‘Evaluating the performance of the lee-carter method for forecasting mortality’, Demography 38(4), 537–549.
  * Li and O’Hare (2017) Li, H. and O’Hare, C. (2017), ‘Semi-parametric extensions of the Cairns-Blake-Dowd model: A one-dimensional kernel smoothing approach’, Insurance: Mathematics and Economics 77, 166 – 176.
  * Li et al. (2015) Li, H., O’Hare, C. and Zhang, X. (2015), ‘A semiparametric panel approach to mortality modeling’, Insurance: Mathematics and Economics 61, 264 – 270.
  * Li and Chan (2005) Li, S.-H. and Chan, W.-S. (2005), ‘Outlier analysis and mortality forecasting: the United Kingdom and Scandinavian countries’, Scandinavian Actuarial Journal 2005(3), 187–211.
  * Lundström and Qvist (2004) Lundström, H. and Qvist, J. (2004), ‘Mortality forecasting and trend shifts: an application of the Lee-Carter model to Swedish mortality data’, International Statistical Review 72(1), 37–50.
  * Park et al. (2009) Park, B. U., Mammen, E., Härdle, W. and Borak, S. (2009), ‘Time series modelling with semiparametric factor dynamics’, Journal of the American Statistical Association 104(485), 284–298.
  * Pena and Box (1987) Pena, D. and Box, G. E. (1987), ‘Identifying a simplifying structure in time series’, Journal of the American Statistical Association 82(399), 836–843.
  * Renshaw and Haberman (2003) Renshaw, A. E. and Haberman, S. (2003), ‘Lee-Carter mortality forecasting with age-specific enhancement’, Insurance: Mathematics and Economics 33(2), 255–272.
  * Stock and Watson (2002) Stock, J. H. and Watson, M. W. (2002), ‘Forecasting using principal components from a large number of predictors’, Journal of the American Statistical Association 97(460), 1167–1179.
  * Stone (1977) Stone, M. (1977), ‘An asymptotic equivalence of choice of model by cross-validation and Akaike’s criterion’, Journal of the Royal Statistical Society. Series B (Methodological) pp. 44–47.
  * Su and Wang (2017) Su, L. and Wang, X. (2017), ‘On time-varying factor models: estimation and testing’, Journal of Econometrics 198(1), 84–101.
  * Tuljapurkar et al. (2000) Tuljapurkar, S., Li, N. and Boe, C. (2000), ‘A universal pattern of mortality decline in the G7 countries’, Nature 405(6788), 789.
  * University of California, Berkeley and Max Planck Institute for Demographic Research (Germany)(n.d.) (USA) University of California, Berkeley (USA) and Max Planck Institute for Demographic Research (Germany) (n.d.), ‘ _Human Mortality Database_ ’.  www.mortality.org [Accessed: 2018.07.10].
  * Yang et al. (2010) Yang, S. S., Yue, J. C. and Huang, H.-C. (2010), ‘Modeling longevity risks using a principal component approach: A comparison with existing stochastic mortality models’, Insurance: Mathematics and Economics 46(1), 254–270.

## Appendix A

Our analysis of the mortality forecasting in the main text focuses on age-
specific US mortality data of the total population. This appendix provides
additional gender-specific results of the estimation and out-of-sample
forecasting for both males and females.

### Male

The results for the male subpopulation are showed in Figure 14 and Figure 15.

In terms of the goodness of fit, the overall MSE for the Lee-Carter model is
$0.008479$, while the MSE for the time-varying model is $0.002679$. From the
perspective of the out-of-sample forecasting precision, the overall MSPE for
the Lee-Carter model is $0.0412585$, while the MSPE for the time-varying model
with the naive method is only $0.02247$, which is approximately $45\%$ less
than the Lee-Carter model. Based on local linear regression, the MSPE of time-
varying model is $0.06222$, which is larger than that for the Lee-Carter
model. Thus, for the male subpopulation, compared with the Lee-Carter model,
the prediction accuracy improves significantly by using the time-varying
factor model based on the naive method. For the short-term forecasting, the
naive method and local regression method have similar performance. However,
for the long-term forecasting, the prediction accuracy of local regression
method deteriorates as time goes by. For the age-specific forecasts, the time-
varying model based on the naive method almost always obtains the best
predictions. And roughly speaking, no matter which method we choose to use,
the time-varying factor model is better at predicting the mortality rates of
the older adulthood ($50\sim 90$).

Figure 14: Year-specific MSPE for both the time-varying model and the Lee-
Carter model; for time-varying model, both naive method and local regression
method are used; Male subpopulation.  Figure 15: Age-specific MSPE for both
the time-varying model and the Lee-Carter model; for time-varying model, both
naive method and local regression method are used; Male subpopulation.

### Female

The results for the female subpopulation are showed in Figure 16 and Figure
17.

In terms of the goodness of fit, the overall MSE for the Lee-Carter model is
$0.007337$, while the MSE for the time-varying model is $0.002061$. From the
perspective of the out-of-sample forecasting precision, the overall MSPE for
the Lee-Carter model is $0.03709$, while the MSPE for the time-varying model
with the naive method is $0.02963$, which is $20\%$ less than the Lee-Carter
model. Based on local linear regression, the MSPE of the time-varying model is
$0.03887$, which is similar to (or slightly larger than) that for the Lee-
Carter model. Thus, for the female subpopulation, compared with the Lee-Carter
model, the prediction accuracy improves significantly by using the time-
varying factor model based on the naive method. For the short-term
forecasting, the naive method and local regression method have similar
performance. However, for the long-term forecasting, the prediction accuracy
of the local regression method deteriorates as time goes by. For the age-
specific forecasts, the time-varying model based on the naive method almost
always obtains the best predictions, except for the age group $40\sim 55$. And
roughly speaking, no matter which method we use, the time-varying factor model
is better at predicting the mortality rates of the young children ($0\sim
10$), young adulthood ($20\sim 40$) and the older adulthood ($55\sim 80$).

Figure 16: Year-specific MSPE for both the time-varying model and the Lee-
Carter model; for time-varying model, both the naive method and local
regression method are used; Female subpopulation.  Figure 17: Age-specific
MSPE for both the time-varying model and the Lee-Carter model; for time-
varying model, both naive method and local regression method are used; Female
subpopulation.

## Appendix B

Figure 18: Plots of the total sum of squared residuals (SSR) versus the length
($k$) of the short-term forecast horizon (based on the hybrid forecasting
method of time-varying factor model); length of forecast horizon: $15$, $20$,
$25$ and $30$

## Appendix C

In Section 5.1, we only present the estimation results of the time-varying
model with a single factor, as the number of factors for the US mortality data
is estimated as 1. This appendix provides some further empirical results of
fitting multi-factor time-varying model to the US mortality data. The number
of factors is 6 (i.e. R = 6) in the following analysis. To solve the
identification problem in the multi-factor context, we require the invertible
matrix ${\bm{H}}_{t}$ as a diagonal matrix and normalize the ${\bm{b}}_{x,t}$
to sum to unity for each $t$. For more details, please refer to Section 2.1.

Figure 19: Plots of the estimated factors for the time-varying factor model
($R=6$)

All six factors are depicted in Figure 19. And in Figure 20, we plot all the
corresponding time-varying factor loadings for ages 20, 40 and 60. We can see
that the first factor and its corresponding factor loadings remain the same as
the results in Section 5.1. The explanation is that theoretically, increasing
the number of factors $R$ does not affect the estimation of the first common
factor $\bm{\kappa}_{t}$ and the corresponding factor loading ${\bm{b}}_{x,t}$
in the localized PCA approach (see the estimation method in Section 2.2). The
first factor captures the major trend of the mortality data, while the
remaining factors are much less important with little information. It
validates our decision to use one factor in the empirical analysis in Section
5.

Figure 20: Plots of the estimated time-varying factor loadings for age 20, 40
and 60 ($R=6$)