\documentclass[preprint,12pt]{elsarticle}
%\documentclass[3p]{elsarticle}

% Natbib setup for author-year style
\usepackage{natbib}
\bibpunct{(}{)}{;}{a}{,}{,}

\usepackage[utf8]{inputenc}
%\VignetteIndexEntry{StMoMo: An R Package for Stochastic Mortality Modelling}
%\VignetteEngine{knitr::knitr}
%!\VignetteEncoding{UTF-8}

%\usepackage[usenames,dvipsnames]{xcolor}
%\usepackage{amssymb}
%\usepackage{subfigure}
%\usepackage{float}
\usepackage[caption = false]{subfig}

%\usepackage{authblk}
%\usepackage{graphicx,epstopdf}
\usepackage{amsmath}
%\usepackage{latexsym}
\usepackage{amsfonts}
%\usepackage{epsfig}
%\usepackage{color}
\usepackage{graphics}
\usepackage{url}
\usepackage[breaklinks=true]{hyperref}
%\usepackage{algorithm, algorithmic, lineno, lscape, ctable}
\usepackage{ctable}
%\usepackage{array}
%\usepackage{rotating}
\usepackage{amsthm}
\usepackage{pdflscape}
\usepackage[margin=1in]{geometry}  %%take this out before submission

% switch off unnnecessary equation numbering
\usepackage{mathtools}
\mathtoolsset{showonlyrefs=true}

\newcommand{\X}{{\cal X}}
\newcommand{\T}{{\cal T}}
\newcommand{\G}{{\cal G}}
\newcommand{\LogL}{{\cal L}}
\newcommand{\Exp}{\mathbb{E}}
\newcommand{\LE}{{\text{\textsl{\r{e}}}}}
\newcommand{\slfrac}[2]{\left.#1\middle/#2\right.}
\DeclareMathOperator{\logit}{logit}

\newcommand{\set}[1]{\left\{ #1 \right\}}
\newcommand{\brac}[1]{\left( #1 \right)}
\renewcommand{\theenumi}{(\alph{enumi})}
\renewcommand{\labelenumi}{\theenumi}


\newcommand{\rcode}[1]{\texttt{#1}}

\newcommand{\rpkg}[1]{\textbf{#1}}
\newcommand{\R}{\rpkg{R} }

\theoremstyle{definition}
\newtheorem{example}{Example}%[subsection]

\newsavebox{\bigleftbox} %Necessary for vertical subfigurees

%Notes in figures or tables
\newcommand\fnote[1]{{\footnotesize{\textit{Note:} {#1}}}}

%Control the position of floats in the document
\setcounter{topnumber}{2}
\setcounter{bottomnumber}{2}
\setcounter{totalnumber}{4}
\renewcommand{\topfraction}{0.85}
\renewcommand{\bottomfraction}{0.85}
\renewcommand{\textfraction}{0.15}
\renewcommand{\floatpagefraction}{0.7}

%Remove Preprint submitted to Elsevier
\makeatletter
\def\ps@pprintTitle{%
  \let\@oddhead\@empty
  \let\@evenhead\@empty
  \def\@oddfoot{\footnotesize\itshape
       \hfill\today}%
  \let\@evenfoot\@oddfoot
}
\makeatother


%==============
% Track changes
%==============
% \newcommand{\draft}[1]{{\color{blue} #1}}
% \newcommand{\note}[1]{{\color{red} #1}}
% \newcommand{\revI}[1]{{\color{Violet} #1}}
% \newcommand{\pietro}[1]{{\color{MidnightBlue} #1}}
% \newcommand{\vlad}[1]{{\color{ForestGreen} #1}}

\newcommand{\draft}[1]{{#1}}
\newcommand{\note}[1]{{#1}}
\newcommand{\revI}[1]{{#1}}
\newcommand{\pietro}[1]{{#1}}
\newcommand{\vlad}[1]{{#1}}


%\newcommand{\revAndrew}[1]{{\color{OrangeRed} #1}}
%\newcommand{\revise}[1]{{\color{Mulberry} #1}}
%\newcommand{\shorten}[1]{{\color{violet} #1}}
%\newcommand{\shorten}[1]{{#1}}

%\newcommand{\draft}[1]{{#1}}
%\newcommand{\revI}[1]{{#1}}
%\newcommand{\revSteve}[1]{{#1}}
%\newcommand{\revAndrew}[1]{{#1}}
%\newcommand{\revise}[1]{{#1}}
%\newcommand{\shorten}[1]{{\color{violet} #1}}
%\newcommand{\shorten}[1]{{#1}}

%\newcommand{\draft}[1]{{#1}}
%\newcommand{\revI}[1]{{#1}}
%\newcommand{\revSteve}[1]{{#1}}
%\newcommand{\revAndrew}[1]{{#1}}
%\newcommand{\revise}[1]{{#1}}
%\newcommand{\shorten}[1]{{\color{violet} #1}}
%\newcommand{\shorten}[1]{{#1}}




%\newcommand{\new}[1]{{#1}}
%\newcommand{\draft}[1]{{#1}}


%\addtolength{\textwidth}{1.5in}
%\addtolength{\hoffset}{-0.75in}
%\addtolength{\textheight}{cin}
%\addtolength{\voffset}{-cin}


%%%%%%%%%%%%%%%%
\begin{document}
%\SweaveOpts{concordance=TRUE}
%%%%%%%%%%%%%%%%

%Hook to control the line width in ouputs
<<lineWidthHook, eval=TRUE,echo = FALSE>>=
library(knitr)
hook_output = knit_hooks$get('output')
knit_hooks$set(output = function(x, options) {
  # this hook is used only when the linewidth option is not NULL
  if (!is.null(n <- options$linewidth)) {
    x = knitr:::split_lines(x)
    # any lines wider than n should be wrapped
    if (any(nchar(x) > n)) x = strwrap(x, width = n)
    x = paste(x, collapse = '\n')
  }
  hook_output(x, options)
})
@


\begin{frontmatter}

% Outcomment only when entries are known. Otherwise leave as is and
%   default values will be used.
%\setcounter{page}{1}
%\VOLUME{00}%
%\NO{0}%
%\MONTH{Xxxxx}% (month or a similar seasonal id)
%\YEAR{0000}% e.g., 2005
%\FIRSTPAGE{000}%
%\LASTPAGE{000}%
%\SHORTYEAR{00}% shortened year (two-digit)
%\ISSUE{0000} %
%\LONGFIRSTPAGE{0001} %
%\DOI{10.1287/xxxx.0000.0000}%

% Enter the full title:
\title{StMoMo: An R Package for Stochastic Mortality Modelling}
\journal{}
%Authors

\author[1]{Andr\'es M. Villegas\corref{cor1}}
\ead{a.villegas@unsw.edu.au}
\author[2]{Pietro Millossovich}
\ead{Pietro.Millossovich.1@city.ac.uk}
\author[2]{Vladimir Kaishev}
\ead{Vladimir.Kaishev.1@city.ac.uk}

\address[1]{ARC Centre of Excellence in Population Ageing Research (CEPAR), UNSW Business School, University of New South Wales, Sydney, Australia}
\address[2]{Cass Business School, City University London, United Kingdom}
\cortext[cor1]{Corresponding author}

\begin{abstract}
{In this paper we mirror the framework of generalised \vlad{(non-)}linear models to define the family of generalised Age-Period-Cohort stochastic mortality models which encompasses the vast majority of stochastic mortality projection models proposed to date\revI{, including the well-known Lee-Carter and Cairns-Blake-Dowd models}. We also introduce the \R package \rpkg{StMoMo} which exploits the unifying framework of the generalised Age-Period-Cohort family to provide tools for fitting stochastic mortality models, \revI{assessing} their goodness of fit and performing mortality projections. We illustrate some of the capabilities of the package by \vlad{performing} a comparison of several stochastic mortality models applied to the England and Wales population.
}

% {Recently there has been a wealth of stochastic mortality models proposed in the demographic and actuarial literature. We introduce the family of generalised Age-Period-Cohort stochastic mortality models which encompasses the vast majority of mortality projection models proposed to date. The \R package \rpkg{StMoMo} exploits the unifying framework of the generalised Age-Period-Cohort family to provide tools for fitting stochastic mortality models, analysing their goodness of fit and performing mortality projections with them. This paper provides an overview of the design and usage of \rpkg{StMoMo}. We illustrate some of the capabilities of the package by doing a comparison of several stochastic mortality models applied to the England and Wales population.
% }

\end{abstract}

% Sample
%\KEYWORDS{Mortality modelling; multipopulation models; socio-economic circumstances, annuity pricing}

% Fill in data. If unknown, outcomment the field
\begin{keyword}
Mortality modelling; mortality forecasting; generalised linear models; \revI{generalised non-linear models}
\end{keyword}
\end{frontmatter}

\section{Introduction}\label{sec:intro}

During the last two centuries developed \pietro{countries} experienced a persistent increase in life expectancy. For instance, \citet{Oeppen2002} estimate that during the last 160 years the world record in female life expectancy at birth has increased at an approximate steady pace of 3 months per year. This increase in life expectancy, though a sign of social progress, poses a challenge to governments, private pension plans and life insurers \revI{because of its} impact on pension and health costs.  Actuaries and demographers have recognised the \revI{problems caused by} an ageing population and rising longevity and have thus devoted \revI{significant attention} to the development of statistical techniques for the modelling and projection of mortality rates.\\  
 
One of the most influential approaches \vlad{to} the stochastic modelling of mortality rates is the parsimonious mortality model proposed by \cite{Lee1992}. This model uses principal component analysis to decompose the \revI{age-time matrix of} mortality rates into a bilinear combination of age and period parameters, with \vlad{the latter} being treated as time series to produce mortality projections. The Lee-Carter model has inspired numerous variants and extensions. For instance, \citet{Lee2001}, \citet{Booth2002}, and  \citet{Brouhns2002} have proposed alternative estimation approaches in order to improve the goodness-of-fit and \vlad{the} forecasting properties of the model.  In particular, \cite{Brouhns2002} propose a more formal statistical approach to estimating the parameters by embedding the Lee-Carter model into a Poisson regression setting. Other authors have extended the Lee-Carter model by including additional terms, such as \revI{multiple} bilinear age-period \pietro{components} \citep{Renshaw2003a, Hyndman2007}, or a cohort effect \revI{term} \citep{Renshaw2006, Currie2006}.\\
  
The two factor \revI{Cairns-Blake-Dowd (CBD)} model introduced by \cite{Cairns2006a} is one of the most prominent \revI{variants of} the Lee-Carter \vlad{model}. \revI{The CBD} model relies on the linearity of the logit of one-year death \pietro{probabilities} at older ages. Specifically, \revI{it} assumes that for a \vlad{given} year the logit of the one-year death probability is a linear function of age, and treats the intercept and slope parameters across years as stochastic processes. \cite{Cairns2009} consider three extensions to the original CBD model by incorporating combinations of a quadratic age term and a cohort effect \vlad{term}. \cite{Plat2009a} has combined \vlad{features of the CBD and the Lee-Carter models} to produce a model that is suitable for full age ranges and captures the cohort effect.\\

Given the abundance and rapid increase in the number of stochastic mortality models proposed in the literature, there have been some recent attempts to find the commonalities among these models. \citet{Hunt2015b} review the structure of mortality models and describe an Age-Period-Cohort model structure which encompasses the vast majority of stochastic mortality models. \citet{Currie2014} shows that many mortality models can be expressed in terms of generalised linear models or generalised non-linear models.\\

In this paper, we build \vlad{upon} the \vlad{works} of \citet{Hunt2015b} and \citet{Currie2014} to define the family of Generalised Age-Period-Cohort stochastic mortality models by mirroring the terminology of generalised linear models. \revI{We also} introduce the \R package \rpkg{StMoMo}\footnote{The acronym StMoMo, pronounced Saint Momo, stands for Stochastic Mortality Modelling. Momo is the king of Carnivals in numerous Latin American festivities \citep{Wikipedia2015}.} which exploits the unifying framework of the Generalised Age-Period-Cohort family to provide computational tools for implementing many of the stochastic mortality models proposed to date. The \rpkg{StMoMo} package is available at \url{http://CRAN.R-project.org/package=StMoMo}. Version 0.3.1 has been used for this paper.\\

Several packages for mortality modelling are available in the \R environment \citep{RCoreTeam2014}. The package \rpkg{demography} \citep{Hyndman2014}, whose usage is explained in detail in \cite{Booth2014}, implements the original Lee-Carter model along with the \cite{Lee2001}, \cite{Booth2002}, and \cite{Hyndman2007} variants. The \rpkg{ilc} package \citep{Butt2014} implements the \citet{Renshaw2006} cohort extension of the Lee-Carter model together with the Lee-Carter model under a Poisson regression framework. The \rpkg{LifeMetrics} \R functions implement the original CBD model and the three extended CBD models considered in \citet{Cairns2009}, along with the Lee-Carter model (using Poisson maximum likelihood), the Age-Period-Cohort model of \citet{Currie2006} and the \citet{Renshaw2006} model. This package, which is not on CRAN, is available at  \url{http://www.macs.hw.ac.uk/~andrewc/lifemetrics/}.\\ 

There are however several \vlad{drawbacks of the existing packages} which our package \rpkg{StMoMo} seeks to overcome. First, the existing packages are based on model-specific fitting algorithms limiting the models available to those already predefined in the packages. By contrast, \rpkg{StMoMo} allow\vlad{s} users to easily expand the number of models available. In addition, whilst \rpkg{StMoMo} provides forecasting and simulation functions for any model within the Generalised Age-Period-Cohort family, the existing packages only provide such \vlad{functions} for a limited number of models. For instance,  \vlad{the} package \rpkg{ilc} only includes forecasting functions for the Lee-Carter model. Similarly, \vlad{simulation with the package \rpkg{LifeMetrics} is limited to the Lee-Carter and the standard CBD models}. Finally, \rpkg{StMoMo} provides \vlad{functions} which are not available in existing packages, such as tools for analysing the goodness-of-fit\footnote{\revI{We note that \rpkg{ilc} also provides some graphical tools for assessing the goodness of fit of the models implemented in that package.}} and evaluating the impact of parameter uncertainty using bootstrapping techniques.\\ 

\pietro{\rpkg{StMoMo} comes with \vlad{a set of} functions for defining an abstract model --- specifying for instance the number of period terms, whether coefficients are parametric or not --- and for fitting a given model. This is particularly useful when estimating several models on a given dataset or a given model to different datasets. The package also provides preset functions for defining the most common models available in the mortality forecasting literature. \vlad{In addition, other models preferred by the user can be created} in a very simple fashion, see Section \ref{sec:GAPC_StMOMO} where several examples are given. Therefore, the flexibility of the package allows a user to quickly build up a battery of different models, and this is particularly useful when seeking the most appropriate mortality model, comparing different models or assessing model risk. \rpkg{StMoMo} is particularly appealing for actuaries managing life and pensions portfolios exposed to longevity risk. The code backing most functions \vlad{implemented} in the package has been extensively used and tested for the development of multi population mortality models for assessing basis risk in longevity risk transactions, see \cite{Haberman2014}.}\\

In the this paper we describe the statistical framework underlying \rpkg{StMoMo} and illustrate \vlad{its usage}. \vlad{For this purpose}, we use as a running example a comparison of several stochastic mortality models fitted to the England and Wales population. This example is in the spirit of the comparison exercises of \cite{Cairns2009,Cairns2011a}, \cite{Haberman2011a} and \cite{Lovasz2011}, \vlad{allowing us to} show how several of the analysis performed in these papers can easily be replicated using \rpkg{StMoMo}. The structure of \vlad{the} paper is as follows. In Section \ref{sec:Notation} we \vlad{introduce} our notation. In Section \ref{sec:GAPC}, we mirror the framework of generalised linear models to define the family of Generalised Age-Period-Cohort (GAPC) stochastic mortality models and \vlad{demonstrate that} many of the mortality models discussed in the literature can be framed within this family. \vlad{In} Section \ref{sec:GAPC_StMOMO} \vlad{we explain how the GAPC family of models is implemented in \rpkg{StMoMo}}. In Section \ref{sec:Fitting}, we describe the fitting of GAPC mortality models and illustrate how this can be accomplished using \rpkg{StMoMo}. \vlad{In Section \ref{sec:GoF} we consider} the evaluation of the goodness-of-fit of GAPC models. In Section \ref{sec:Forecasting} we discuss the forecasting and simulation of GAPC models using time series techniques. Section \ref{sec:ParameterUncertainty} describes the use of bootstrapping techniques to incorporate parameter uncertainty in the estimation and forecasting of GAPC mortality models. Finally, \vlad{in Section \ref{sec:Conclusions} we provide some conclusions and discuss  possible extensions of the \rpkg{StMoMo} package}. 




%   \begin{itemize}  
% %  \item \draft{Review of mortality modelling}
% %   \begin{itemize} 
% %     \item \draft{Review of papers comparing models. \citep{Haberman2011a, Cairns2009, Lovasz2011, VanBerkum2014}}
% %    \end{itemize}
% %   \item \draft{Review of mortality modelling in R}	
% %   \begin{itemize}
% %     \item \draft{Demography}
% %     \item \draft{ilc \cite{Butt2014}}
% %     \item \draft{Lifemetrics}	
% %     \item \draft{MortaltySmooth \cite{Camarda2012}}
% %   	\item \draft{ELT}
% % 		\item \draft{Paper by Currie \cite{Currie2014}}
% %   \end{itemize}
% %   
% %   \item \draft{Our contributions}  
% %   \begin{itemize}
% %     \item \draft{Unifying framework for mortality modelling (GLM formulations)}
% %     \item \draft{Exploit this unifying framework to produce a general package for mortality modelling.}
% %   \end{itemize}
%   
%   \item \draft{Justification of the package and what differes from existing R packages.}
%   \begin{itemize}
%   \item \draft{Flexibility of running different models under the same environment (unifying framework)}
% 	\item \draft{Provide tools for comparing models and selecting the most appropriate}
% 	\item \draft{Cover model families which have not been implemented in previously existing packages}
% 	\item \draft{Possibility extending the package by implementing new models under the same framework as the literature progress} 
% 	\item \draft{First step in developing multi-population mortality models (see Hymans report)}
% 	\item \draft{Useful in the lab for teaching students}
%   \item \draft{Fitting based on \rpkg{gnm}. This gives the user flexibility to control parameters in the fitting algorithm (e.g. tolerance, intial parameter values, etc.)}
%  \end{itemize}
% %\item \draft{Objective of the paper: to introduce the package its main functions and to illustrate its usage on various examples.}
%    
% \end{itemize}



\section{Notation and Data}\label{sec:Notation}

Let the random variable $D_{xt}$ denote the number of deaths in a population at age $x$ last birthday during calendar year $t$. Also let $d_{xt}$ denote the observed number of deaths, $E^c_{xt}$ the central exposed to risk at age $x$ in year $t$, and $E^0_{xt}$ the corresponding initial exposed to risk. The one-year mortality rate for an individual aged $x$ last birthday and in calendar year $t$, denoted $q_{xt}$, can be estimated as $\hat{q}_{xt}=\slfrac{d_{xt}}{E^0_{xt}}$. The force of mortality and central death rates are denoted by $\mu_{xt}$ and $m_{xt}$, respectively, with the empirical estimate of the latter being $\hat{m}_{xt}=\slfrac{d_{xt}}{E^c_{xt}}$. Under the assumption that the force of mortality is constant over each year of age \pietro{and calendar year}, i.e. from age $x$ to age $x+1$ and $t$ to $t+1$, then the force of mortality $\mu_{xt}$ and the death rate $m_{xt}$ coincide. We assume that this is the case throughout.\\

%Let $D_{xt}$ denote the number of death age $x$ last birthday calendar during calendar year $t$, $E^c_{xt}$ the central exposed to risk at age $x$ in year $t$ and $E^0_{xt}$ the corresponding initial exposed to risk at age $x$ in year $t$. The 1-year mortality rate for an individual aged $x$ last birthday and in calendar year $t$, denoted $q_{xt}$, can be computed as $q_{xt}=\slfrac{D_{xt}}{E^0_{xt}}.$ The corresponding central death rate, $m_{xt}$, is defined as $m_{xt}=\slfrac{D_{xt}}{E^c_{xt}}$. The force of mortality $\mu_{xt}$ is the instantaneous death rate at exact time $t$ and exact age $x$. We note that under the assumption that the force of mortality is constant over each year of age, i.e. from age $x$ to age $x+1$, then the force of mortality $\mu_{xt}$ and the death rate $m_{xt}$ coincide.\\

In \rpkg{StMoMo} and throughout this paper we assume that deaths, $d_{xt}$, and either central exposures, $E^c_{xt}$, or initial exposures, $E^0_{xt}$, are available in a rectangular array format comprising ages (on the rows) $x = x_1, x_2, \ldots, x_k$, and calendar years (on the columns) $t = t_1,t_2, \ldots, t_n$, When only central exposures are available and initial exposures are required (or vice-versa), one can approximate the initial exposures by adding half the matching reported numbers of deaths to the central exposures, i.e, $E^0_{xt}\approx E^c_{xt} + \frac{1}{2}d_{xt}$. When the context is clear, we may write $E_{xt}$ to refer to $E^0_{xt}$ or $E^c_{xt}$.

\section{Generalised Age-Period-Cohort stochastic mortality models}\label{sec:GAPC}
Some authors have recently sought to identify the similarities amongst stochastic mortality models. For instance, \citet{Hunt2015b} describe an Age-Period-Cohort model structure which encompasses the vast majority of stochastic mortality models. In another interesting contribution, \citet{Currie2014} shows \vlad{that} many common mortality models can be expressed in the standard terminology of generalised linear \pietro{or non-linear models}. In this section, we build \vlad{upon} the aforementioned papers to define the family of Generalised Age-Period-Cohort (GAPC) stochastic mortality models. \\

Akin to generalised linear models (see e.g. \citet{McCullagh1989}), a GAPC stochastic mortality model is comprised of four components:

\begin{enumerate}[i.]
 %Random Component
 \item The \textit{random component}: \pietro{the numbers of deaths $D_{xt}$} \revI{follow a Poisson distribution or a Binomial distribution}\footnote{\revI{More precisely, conditionally on $(\mu_{xt})$, the numbers of deaths are independent and follow a Poisson distribution. Similarly, in the Binomial case, the numbers of deaths are independent conditionally on $(q_{xt})$.}}, so that 
 \begin{equation*}\label{eq:Poisson}
D_{xt}\sim \mathrm{Poisson}(E^c_{xt}\mu_{xt})
\end{equation*}
or 
\begin{equation*}\label{eq:Binomical}
D_{xt}\sim \mathrm{Binomial}(E^0_{xt}, q_{xt}),
\end{equation*}
with $\Exp\left(\slfrac{D_{xt}}{E^c_{xt}}\right) = \mu_{xt}$ and $\Exp\left(\slfrac{D_{xt}}{E^0_{xt}}\right) = q_{xt}$, respectively.\\ 
 
 
 %Systematic Component
 \item The \textit{systematic component}: following \citet{Hunt2015b} the effects of age $x$, calendar year $t$  and year-of-birth (cohort) $c=t-x$ are captured through a \textit{predictor} $\eta_{xt}$ given by:
 \begin{equation*}\label{eq:predictor}
  \eta_{xt} = \alpha_x + \sum_{i=1}^N \beta_x^{(i)}\kappa_t^{(i)} + \beta_x^{(0)}\gamma_{t-x}.
 \end{equation*}
Here:
\begin{itemize}
  \item The term $\alpha_x$ is a static age function capturing the general shape of mortality by age.
  \item $N\geq 0$ is an integer indicating the number of age-period terms describing the mortality trends, with each time index $\kappa_t^{(i)}$, $i = 1,\ldots, N$, contributing in specifying the mortality trend and $\beta_x^{(i)}$ modulating its effect across ages.
 \item The term $\gamma_{t-x}$ accounts for the cohort effect  with $\beta_x^{(0)}$ modulating its effect across ages.
\end{itemize}
The age modulating terms $\beta_x^{(i)}$, $i = 0,1,\ldots, N$, can be either pre-specified functions of age, i.e. $\beta_x^{(i)}\equiv f^i(x)$, as in CBD type models, or non-parametric terms without any prior structure and which need to be estimated as in the Lee-Carter model. In the GAPC family we assume that the period indexes $\kappa_t^{(i)}$, $i = 1,\ldots, N$, and the cohort index $\gamma_{t-x}$ are stochastic processes rather than smooth functions of time or cohort. This is the key feature that allows the stochastic projection of GAPC models and thus the generation of probabilistic forecasts of future mortality rates.\\

%Link function 
\item The \textit{link function} \revI{$g$} associating the random component and the systematic component so that 
 \begin{equation*}\label{eq:link}
  g\left(\Exp\left(\frac{D_{xt}}{E_{xt}}\right)\right) = \eta_{xt}  .
 \end{equation*}
Although a number of link functions would be possible it is convenient to use the so-called canonical link and pair the Poisson distribution with the log link function and the Binomial distribution with the logit link function (see e.g. \citet{Currie2014} for a discussion of this in the context of mortality models and \citet{McCullagh1989} in the wider context of GLMs).\\

%constraint function 
\item The \textit{set of parameter constraints}: most stochastic mortality models are only identifiable up to a transformation and thus require parameter constraints to ensure unique parameter estimates. These parameter constraints are applied \pietro{through a} \textit{constraint function} \revI{$v$} which maps an arbitrary \revI{vector} of \vlad{parameters}
\[\theta := \left(\alpha_{x}, \beta_x^{(1)},..., \beta_x^{(N)}, \kappa_t^{(1)},..., \kappa_t^{(N)}, \beta_x^{(0)}, \gamma_{t-x}\right)\] 
into a \revI{vector} of \vlad{transformed parameters} \revI{\[v(\theta)=\tilde{\theta} = \left(\tilde{\alpha}_{x}, \tilde{\beta}_x^{(1)},..., \tilde{\beta}_x^{(N)}, \tilde{\kappa}_t^{(1)},..., \tilde{\kappa}_t^{(N)}, \tilde{\beta}_x^{(0)}, \tilde{\gamma}_{t-x}\right)\]}
satisfying the model constraints \revI{with no effect on the predictor $\eta_{xt}$ (i.e. $\theta$ and $\tilde{\theta}$ result in the same $\eta_{xt}$)}.
\end{enumerate}

Most stochastic mortality models proposed in the literature belong to the GAPC family. This includes the original Lee-Carter model, the extensions of the Lee-Carter proposed in \citet{Renshaw2003a, Renshaw2006}, the original CBD model, and the extended CBD models of \cite{Cairns2009}. In addition, all the model structures considered in \citet{Haberman2011a}, \citet{Lovasz2011}  and \citet{VanBerkum2014}, as well as the models of \cite{Plat2009a}, \cite{Aro2011}, \cite{OHare2012}, \citet{Borger2013} and \citet{Alai2014b}, are part of the GAPC class of models.\footnote{\revI{We note however that models which rely on the smoothness of mortality over both age and time, such as the graduation approach of \citet{Renshaw1996} and the P-Spline model of \citet{Currie2006}, do not belong to the GAPC family. The P-Spline approach is implemented in the \rpkg{MortalitySmooth} \citep{Camarda2012} \R package.}}\\

Next, we describe in detail some of these models highlighting how they can  be framed within the GAPC family.

%\draft{There are however some models which do not lie within ghe GAPC family. P-Splines \citep{Currie2004, Camarda2012} and \citet{Sithole2000}}.\\ 

% Most stochastic mortality models proposed in the literature belong to the GAPC family. This includes all the Lee-Carter type models and CBD type models discussed in sections \ref{sec:LeeCarterReview} and \ref{sec:CBDReview}, but also some other models. In particular, all the model structures considered in \citet{Haberman2011a} and \citet{VanBerkum2014}, as well as the models of \cite{Aro2011} and \citet{Borger2013}, are part of the GAPC family.\\

% Next, we show how sevaral common stochastic mortality models can be framed within the GAPC family.

\subsection{Lee-Carter model under a Poisson setting}\label{ex:LC}
%\begin{example}[Lee-Carter model under a Poisson setting]\label{ex:LC}
The Lee Carter model as implemented by \citet{Brouhns2002} assumes a Poisson distribution of the deaths using a log \textit{link function} to target the force of mortality $\mu_{xt}$. The predictor structure proposed by \citet{Lee1992} assumes that there is a static age function, $\alpha_x$, a unique non-parametric age-period term ($N=1$), and no cohort effect. Thus, the predictor is given by:
\begin{equation}\label{eq:LCpredictor}
\eta_{xt} = \alpha_x + \beta^{(1)}_x \kappa^{(1)}_t
\end{equation}
\vlad{In order to project mortality, the time index $\kappa^{(1)}_t$ is modelled and forecasted using ARIMA processes. Typically, a random walk with drift {has been shown to provide} a reasonable fit, that is,
\begin{equation}\label{eq:LCkappa}
\kappa^{(1)}_t = \delta + \kappa^{(1)}_{t-1} + \xi_t, \qquad \xi_t\sim N(0,\sigma_{\kappa}), 
\end{equation}
where $\delta$ is the drift parameter and $\xi_t$ is Gaussian white noise process with variance $\sigma_{\kappa}$.}

The Lee-Carter model is only identifiable up to a transformation, as for arbitrary real constants $c_1$ and $c_2\neq 0$ the parameters in Equation \eqref{eq:LCpredictor} can be transformed in the following way 
\begin{equation}\label{eq:LCTransformations}
\left(\alpha_x,\beta^{(1)}_x, \kappa^{(1)}_t\right)\to\left(\alpha_x + c_1\beta^{(1)}_x, \frac{1}{c_2}\beta^{(1)}_x, c_2(\kappa^{(1)}_t - c_1)\right),
\end{equation}
leaving $\eta_{xt}$ unchanged. To ensure identifiability of the model \citet{Lee1992} suggest the  following \textit{set of parameter constraints}
\begin{equation}\label{eq:LCConstraints}
\sum_x \beta^{(1)}_x = 1, \qquad \sum_t \kappa^{(1)}_t = 0, 
\end{equation}
which can be imposed by choosing
\begin{equation}\label{eq:LCConstants}
c_1 = \frac{1}{n}\sum_t \kappa^{(1)}_t, \qquad c_2 = \sum_x \beta^{(1)}_x
\end{equation}
in transformation \eqref{eq:LCTransformations}.

%\end{example}

\subsection{Renshaw and Haberman model: Lee-Carter with cohort effects}\label{ex:RH}
\cite{Renshaw2006} generalise the Lee-Carter model by incorporating a cohort effect to obtain the predictor:
\begin{equation}\label{eq:LCCohortEffect}
\eta_{xt} = \alpha_x + \beta^{(1)}_x \kappa^{(1)}_t   + \beta^{(0)}_x \gamma_{t-x}
\end{equation} 
\vlad{Mortality projections for this model are derived using time series forecast of the estimated $\kappa^{(1)}_t$ and $\gamma_{t-x}$, generated using univariate ARIMA processes under the assumption of independence between the period and the cohort effects.}\\

To estimate the model \cite{Renshaw2006} assume a Poisson distribution of deaths (\textit{random component}) and use a log \textit{link} function targeting the force of mortality $\mu_{xt}$. As with the Lee-Carter model, the predictor $\eta_{xt}$ is invariant to the transformation:
\begin{align}
\left(\alpha_x,\beta^{(1)}_x, \kappa^{(1)}_t, \beta^{(0)}_x, \gamma_{t-x}\right)\to &\bigg(\alpha_x + c_1\beta^{(1)}_x + c_2\beta^{(1)}_x , \frac{1}{c_3}\beta^{(1)}_x, \notag\\
& c_3(\kappa^{(1)}_t - c_1), \frac{1}{c_4}\beta^{(0)}_x,  c_4(\gamma_{t-x} - c_2)\bigg),\label{eq:RHtransformations}
\end{align}
where $c_1$, $c_2$, $c_3\neq 0$ and $c_4\neq 0$ are real constants. Identifiability of the model can be ensured using the following \textit{set of parameter constraints}:
\begin{equation}\label{eq:RHGAPCconstraints}
\sum_x \beta^{(1)}_x = 1, \qquad \sum_t \kappa^{(1)}_t = 0, \qquad \sum_x \beta^{(0)}_x = 1, \qquad \sum_{c=t_1-x_k}^{t_n - x_1}\gamma_c = 0, 
\end{equation}
which can be imposed by setting
\begin{equation}\label{eq:RHGAPCconsFun}
c_1 = \frac{1}{n}\sum_t \kappa^{(1)}_t, \quad c_2 = \frac{1}{n+k-1}\sum_{c=t_1-x_k}^{t_n - x_1}\gamma_c, \quad c_i = \sum_x \beta^{(i)}_x, \quad i = 3, 4,
\end{equation}
in transformation \eqref{eq:RHtransformations}.\\

\cite{Renshaw2006} also consider several substructures of the predictor \eqref{eq:LCCohortEffect} obtained by setting to a constant one or both of the age modulating terms. Of particular interest is the substructure \pietro{obtained by setting $\beta_{x}^{(0)} = 1$,}
\begin{equation}\label{eq:CohortEffectH1}
\eta_{xt} = \alpha_x + \beta^{(1)}_x \kappa^{(1)}_t   + \gamma_{t-x},
\end{equation} 
which has been suggested by \cite{Haberman2011a} as a simpler structure that resolves some stability issues of the original model.

\subsection{Age-Period-Cohort model}\label{ex:APC}
Another commonly used substructure of the Renshaw and Haberman model is the so-called Age-Period-Cohort (APC) model, \pietro{corresponding to $\beta_{x}^{(1)} = 1$, $\beta_{x}^{(0)} = 1$,}
\begin{equation}\label{eq:CohortEffectAPC}
\eta_{xt} = \alpha_x +  \kappa^{(1)}_t   + \gamma_{t-x},
\end{equation} 
which has a long-standing tradition in the fields of medicine and demography (see, e.g., \cite{Clayton1987} and \cite{Hobcraft1982}), but has not been widely used in the actuarial literature until \pietro{it was} considered by \cite{Currie2006}. The APC model is known to be \pietro{invariant with respect to} the following two transformations:

\begin{align}
\left(\alpha_x, \kappa^{(1)}_t, \gamma_{t-x}\right) &\to \left(\alpha_x + \phi_1 - \phi_2x, \kappa^{(1)}_t + \phi_2t, \gamma_{t-x} - \phi_1 - \phi_2(t-x)\right)\label{eq:APCTransGamma}\\
\left(\alpha_x, \kappa^{(1)}_t, \gamma_{t-x}\right) &\to \left(\alpha_x + c_1, \kappa^{(1)}_t - c_1, \gamma_{t-x}\right)\label{eq:APCTransKappa},
\end{align}
where $c_1$, $\phi_1$, and $\phi_2$ are real constants. However, we can ensure identifiability of the model by imposing the \textit{set of constraints}:
\begin{equation}\label{eq:APCconstraints}
\sum_t \kappa^{(1)}_t = 0, \qquad \sum_{c=t_1-x_k}^{t_n - x_1}\gamma_c = 0, \qquad \sum_{c=t_1-x_k}^{t_n - x_1}c\gamma_c = 0,
\end{equation}
where the last two constraints \pietro{imply} that the cohort effect fluctuates around zero with no discernible linear trend. Following \citet[Appendix A]{Haberman2011a}, the constraints on the cohort effect can be imposed by applying transformation \eqref{eq:APCTransGamma} with constants $\phi_1$ and $\phi_2$ obtained by regressing $\gamma_{t-x}$ on $t-x$, so that 
\begin{equation}\label{eq:APCConstFun1}
\gamma_{t-x} = \phi_1 + \phi_2(t-x) + \epsilon_{t-x}, \quad \epsilon_{t-x} \sim N(0,\sigma^2) \quad \text{i.i.d.}
\end{equation}
The constraint on the period index can then be imposed by applying transformation \eqref{eq:APCTransKappa} with 
\begin{equation}\label{eq:APCconsFun2}
c_1 = \frac{1}{n}\sum_t \kappa^{(1)}_t.
\end{equation}


\subsection{CBD model}\label{ex:CBD}
%\begin{example}[CBD model under a Binomial setting]
\citet{Cairns2006a} propose a predictor structure with two age-period terms ($N=2$) with pre-specified age-modulating parameters $\beta^{(1)}_x = 1$ and $\beta^{(2)}_x = x-\bar{x}$, no static age function and no cohort effect. Thus, the predictor of the CBD model is given by: 
\begin{equation}\label{eq:CBDGAPCpredictor}
\eta_{xt} =  \kappa_t^{(1)} + (x-\bar{x})\kappa_t^{(2)},
\end{equation}
where $\bar{x}$ is the average age in the data. \vlad{\citet{Cairns2006a} obtain mortality forecasts by projecting the period effects $\kappa_t^{(1)}$ and $\kappa_t^{(2)}$ using a bivariate random walk with drift}.\\

The CBD model does not have identifiability issues and hence the \textit{set of parameter constraints} is empty. In order to estimate the parameter of the CBD model we can follow \citet{Haberman2011a} and  assume a Binomial distribution of deaths using a logit \textit{link function} targeting the one-year death probabilities $q_{xt}$.

\subsection{M7: Quadratic CBD model with cohort effects}\label{ex:M7}
\cite{Cairns2009} extend the original CBD model by adding a cohort effect and a quadratic age effect to obtain the predictor:
\begin{equation}\label{eq:M7}
\eta_{xt} = \kappa_t^{(1)} + (x-\bar{x})\kappa_t^{(2)} + \kappa^{(3)}_t \left((x-\bar{x})^2-\hat{\sigma}^2_x\right) + \gamma_{t-x}, 
\end{equation}
where $\hat{\sigma}^2_x$ is the average value of $(x-\bar{x})^2$. This model, usually referred to as model M7, is not identifiable \pietro{as} the transformation
\begin{align}
\left(\kappa^{(1)}_t, \kappa^{(3)}_t, \kappa^{(3)}_t, \gamma_{t-x}\right) \to& \bigg(\kappa^{(1)}_t + \phi_1 + \phi_2 (t-\bar{x}) + \phi_3\left((t-\bar{x})^2+\hat{\sigma}^2\right), \kappa^{(2)}_t - \phi_2-2\phi_3(t-\bar{x}),\notag \\
 & \kappa^{(3)}_t + \phi_3, \gamma_{t-x} - \phi_1 - \phi_2(t-x) - \phi_3(t-x)^2\bigg),\label{eq:M7Trans}
\end{align}
for real constants $\phi_1$, $\phi_2$ and $\phi_3$ \pietro{, leave the predictor unchanged}. To identify the model \cite{Cairns2009} impose the \textit{set of constraints}:
\begin{equation}\label{eq:M7constraints}
\sum_{c=t_1-x_k}^{t_n - x_1}\gamma_c = 0, \qquad \sum_{c=t_1-x_k}^{t_n - x_1}c\gamma_c = 0, \qquad \sum_{c=t_1-x_k}^{t_n - x_1}c^2\gamma_c = 0, 
\end{equation}
which  ensure that the cohort effect fluctuates around zero and has no discernible linear or quadratic trend. Following \citet[Appendix A]{Haberman2011a}, these constraints can be imposed by applying transformation \eqref{eq:M7Trans} with constants $\phi_1$, $\phi_2$ and $\phi_3$ obtained by regressing $\gamma_{t-x}$ on $t-x$ and $(t-x)^2$, so that 
\begin{equation}\label{eq:M7ConstFun1}
\gamma_{t-x} = \phi_1 + \phi_2(t-x) + \phi_3(t-x)^2 + \epsilon_{t-x}, \quad \epsilon_{t-x} \sim N(0,\sigma^2) \quad \text{i.i.d.}
\end{equation}
\cite{Cairns2009} also consider the simpler predictor structures
\begin{align}
\eta_{xt} &= \kappa_t^{(1)} + (x-\bar{x})\kappa_t^{(2)} + \gamma_{t-x}, \label{eq:M6}\\
\eta_{xt} &= \kappa_t^{(1)} + (x-\bar{x})\kappa_t^{(2)} + (x_c - x)\gamma_{t-x}, \label{eq:M8}
\end{align}
where $x_c$ is a constant parameter to be estimated. These structures are typically referred to as models M6 and M8, respectively. 

\subsection{Plat model}\label{ex:Plat}
%\begin{example}[Plat model for older ages]\label{ex:Plat}
\citet{Plat2009a} combines the CBD model with some features of the Lee-Carter model to produce a model that is suitable for full age ranges and captures the cohort effect. The proposed predictor structure assumes that there is a static age function, $\alpha_x$, three age-period terms ($N=3$) with pre-specified age-modulating parameters $\beta^{(1)}_x = 1$, $\beta^{(2)}_x = \bar{x} -x$, $\beta^{(3)}_x = (\bar{x}-x)^+=\max(0, \bar{x}-x)$, and a cohort effect with pre-specified age-modulating parameters $\beta^{(0)}_x = 1$. Thus, the predictor is given by:  
\begin{equation}\label{eq:Plat}
\eta_{xt} =  \alpha_x  + \kappa_t^{(1)} + (\bar{x}-x)\kappa_t^{(2)} + (\bar{x}-x)^+\kappa^{(3)}_t + \gamma_{t-x}.
\end{equation}
\citet{Plat2009a} targets the force of mortality $\mu_{xt}$ with the log \textit{link} and estimates the parameters of the model by assuming a Poisson distribution of the deaths.  The following parameter transformations leave the predictor in \eqref{eq:Plat} unchanged:
\begin{align}
\left( \alpha_x, \kappa^{(1)}_t, \kappa^{(2)}_t, \kappa^{(3)}_t,  \tilde{\gamma}_{t-x} \right) &\to \Big(\alpha_x + \phi_1 - \phi_2 x + \phi_3 x^2 ,\kappa^{(1)}_t + \phi_2t + \phi_3(t^2-2\bar{x}t), \notag\\
& \kappa^{(2)}_t + 2\phi_3t, \kappa^{(3)}_t, \gamma_{t-x} - \phi_1 - \phi_2(t-x) - \phi_3(t-x)^2\Big) \label{eq:PlatTrans2}\\
\left( \alpha_x, \kappa^{(1)}_t, \kappa^{(2)}_t, \kappa^{(3)}_t, \gamma_{t-x} \right) &\to \Big(\alpha_x + c_1 + c_2 (\bar{x}-x) + c_3(\bar{x}-x)^+,  \notag\\
& \kappa^{(1)}_t - c_1, \kappa^{(2)}_t - c_2, \kappa^{(3)}_t - c_3 , \gamma_{t-x} \Big) \label{eq:PlatTrans1},
\end{align}
where $c_1$, $c_2$, $c_3$, $\phi_1$, $\phi_2$, and $\phi_3$ are any real constants. The following \textit{set of parameter constraints} can be imposed to ensure identifiability:
\begin{equation} \label{eq:PlatConst}
\sum_t \kappa^{(1)}_t = 0,  \sum_t \kappa^{(2)}_t = 0, \sum_t \kappa^{(3)}_t = 0, \sum_{c=t_1-x_k}^{t_n - x_1}\gamma_c = 0, \sum_{c=t_1-x_k}^{t_n - x_1} c\gamma_c = 0, \sum_{c=t_1-x_k}^{t_n - x_1} c^2\gamma_c = 0  
\end{equation}
The first three constraints ensure that the period indexes are centred around zero, while the last three constraints ensure that the cohort effect fluctuates around zero and has no linear or quadratic trend. Following \citet[Appendix A]{Haberman2011a}, the constraints on the cohort effect can be imposed by applying transformation \eqref{eq:PlatTrans2} with constants $\phi_1$, $\phi_2$, and $\phi_3$ obtained by regressing $\gamma_{t-x}$ on $t-x$ and $(t-x)^2$, so that 
\begin{equation}\label{eq:PlatConstFun1}
\gamma_{t-x} = \phi_1 + \phi_2(t-x) + \phi_3(t-x)^2 + \epsilon_{t-x}, \quad \epsilon_{t-x} \sim N(0,\sigma^2) \quad \text{i.i.d. .}
\end{equation}
The constraints on the period indexes can then be imposed by applying transformation \eqref{eq:PlatTrans1} with
\begin{equation}\label{eq:PlatConstFun3}
c_i = \frac{1}{n}\sum_t \kappa^{(i)}_t, \qquad i=1,2,3.
\end{equation}
In the cases where only older ages are of interest, \cite{Plat2009a} suggests to drop the third period term from predictor \eqref{eq:Plat}\footnote{Note that the reduced Plat model is essentially the M6 model of \cite{Cairns2009} with an added static age term $\alpha_x$.}:   
\begin{equation}\label{eq:PlatOld}
\eta_{xt} = \alpha_x + \kappa^{(1)}_t + (\bar{x}-x)\kappa^{(2)}_t + \gamma_{t-x}.
\end{equation}
We note that this reduced Plat model has the same identifiability issues as the complete model with the omission of the transformations and constraints involving $\kappa^{(3)}_t$ and $c_3$. 


%  \subsection{The Lee-Carter model}\label{sec:LeeCarterReview}
%  
% The model of \cite{Lee1992} is one of the most influential models for mortality projection in recent times. Indeed, most of the recent developments in mortality forecasting  build on its basic ideas. \cite{Lee1992} assume that the central death rate $m_{xt}$ follows a log-bilinear form:
% \begin{equation}\label{eq:LC}
% \log m_{xt}= \alpha_x + \beta_x\kappa_t + \epsilon_{xt}                                                                        
% \end{equation}
% where $\alpha_x$ captures the general age-specific mortality pattern, $\kappa_t$ is a time varying mortality index representing the overall level of mortality in year $t$, and $\beta_x$ measures the age-specific response to changes in the general level $\kappa_t$. The error term $\epsilon_{xt}$, with mean $0$ and variance $\sigma^2_{e}$, captures the variations not explained by the model.\\ 
%  
% The calibration of the Lee-Carter model poses two problems. First, the parameters of the model cannot be estimated using simple regression methods as all the variables on the right-hand side of~\eqref{eq:LC} are non-observable. Second, the parameters are only identifiable up to a transformation. Specifically, for arbitrary real constants $c_1$ and $c_2\neq 0$ the parameters in Equation \eqref{eq:LC} can be transformed in the following way 
% \begin{equation}\label{eq:LCtransformations}
% \left\{ \tilde{\alpha}_x, \tilde{\beta}_x, \tilde{\kappa}_t \right\} = \left\{\alpha_x + c_1\beta_x, \frac{1}{c_2}\beta_x, c_2(\kappa_t-c_1) \right\} ,
% \end{equation}
% leaving the fitted mortality rates unchanged. To tackle these problems, \cite{Lee1992} obtain a {least-squares} solution for $\alpha_x$, $\beta_x$ and $\kappa_t$ by means of singular value decomposition techniques and impose the constraints 
% \begin{align}
% \sum_x \beta_x &= 1 \label{eq:ScaleBeta} \\
% \sum_t \kappa_t &= 0 \label{eq:LevelKappa}
% \end{align}
% ensuring the identifiability of the model.\\
% 
% {Instead of keeping $\kappa_t$'s obtained from singular value decomposition, \cite{Lee1992} suggest a second stage process whereby the $\kappa_t$'s are re-estimated so that the fitted number of deaths matches the observed number of deaths in each year. This adjustment gives higher weight to mortality rates at older ages, avoiding the discrepancies induced by modelling mortality rates on a logarithmic scale.} In order to project mortality, the {adjusted} time index $\kappa_t$ is modelled and forecasted using Box-Jenkins time series methods. Typically, a random walk with drift {has been shown to provide} a reasonable fit, that is,
% \begin{equation}\label{eq:LCkappa}
% \kappa_t = d + \kappa_{t-1} + \xi_t,
% \end{equation}
% where $d$ is the drift parameter and $\xi_t$ is an error term.\\
% 
% The Lee-Carter model has inspired numerous variants and extensions. For instance, as opposed to matching the total number of deaths as originally suggested by \cite{Lee1992}, \cite{Lee2001} have proposed to adjust the $\kappa_t$'s to reproduce period life expectancy in each year, while  \cite{Booth2002} have proposed to adjust the $\kappa_t$'s to fit the age distribution of deaths. \cite{Hyndman2007} have extended the Lee-Carter model by adopting a functional data approach. \R implementations of the original Lee-Carter model along with the \cite{Lee2001}, \cite{Booth2002}, and \cite{Hyndman2007} variants are provided by package \rpkg{demography} whose usage is explained in detail in \cite{Booth2014}, \\ 
% 
% \cite{Brouhns2002} point out that the {least-squares} parameter estimation procedure proposed by \cite{Lee1992} assumes that the error terms are homoskedastic. This is an unrealistic assumption as the logarithm of the observed force of mortality tends to be more volatile at higher ages than at younger ages. To deal with this shortcoming \cite{Brouhns2002} embed the Lee-Carter in a Poisson regression setting targeting the force of mortality $\mu_{xt}$. The Poisson approach assumes that  
% 
% \begin{equation}\label{eq:LClogbilinear}
% D_{xt}\sim \mathrm{Poisson}(E^c_{xt}\mu_{xt})
% \end{equation}
% with 
% \begin{equation}\label{eq:LClogbilinearMu}
% \log \mu_{xt} = \alpha_x + \beta_x \kappa_t
% \end{equation}
% where the parameters are subject to the same constraints as in the traditional Lee-Carter model and have the same interpretation. Parameter estimates of the model are obtained by maximising the log-likelihood of the model.\\
% 
% The Lee-Carter model uses exclusively the first principal component of the matrix of centred log-death rates \citep{Bell1997}, resulting in a perfect correlation of mortality improvements across ages. For the purpose of allowing a richer correlation structure in the rates of improvement across ages, \cite{Renshaw2003a} have investigated the feasibility of extending the Lee-Carter model by incorporating the second principal component. Thus, \cite{Renshaw2003a} have proposed the model 
% \begin{equation}\label{eq:LCAgeEnhanced}
%   \log \mu_{xt} = \alpha_x + \beta^{(1)}_x \kappa^{(1)}_t + \beta^{(2)}_x \kappa^{(2)}_t	
% \end{equation}
% where $\kappa^{(1)}_t$ and $\kappa^{(1)}_t$ are period effects. Mortality forecasts are then obtained by projecting the period effects either as separate univariate ARIMA processes \citep{Renshaw2003a} or as bivariate time series \citep{Renshaw2005}. \\
% 
% 
% In some countries it has been observed that mortality improvements depend not only on age and period but on year of birth. In the particular case of the United Kingdom the generations born between 1925 and 1945 have experienced a more rapid mortality improvement relative to those born before or after that period \citep{Willets2004}. Some plausible explanations for this so-called cohort effect include differential dietary habits in early life, the introduction of the Welfare State during the late 1940's, and differential smoking patterns across cohorts. \cite{Renshaw2006} have extended the Lee-Carter model to incorporate the cohort effect. The Renshaw and Haberman extension of the Lee-Carter model can be written as:
% \begin{equation}\label{eq:LCCohortEffect}
% \log \mu_{xt} = \alpha_x + \beta^{(1)}_x \kappa_t   + \beta^{(0)}_x \gamma_{t-x}
% \end{equation} 
% where $\alpha_x$ gives the average age-pattern of mortality, $\kappa_t$ represent the general level of mortality in year $t$, $\gamma_{t-x}$ represents the cohort effect, and $\beta^{(1)}_x$ and $\beta^{(0)}_x$ modulate the age-specific response to changes in $\kappa_t$ and $\gamma_{t-x}$.  Other commonly used substructures of the Renshaw and Haberman model are 
% \begin{equation}\label{eq:CohortEffectH1}
% \log \mu_{xt} = \alpha_x + \beta^{(1)}_x \kappa_t	 + \gamma_{t-x} \quad (\beta_{x}^{(0)} = 1)
% \end{equation} 
% and  
% \begin{equation}\label{eq:CohortEffectAPC}
% \log \mu_{xt} = \alpha_x +  \kappa_t	 + \gamma_{t-x} \quad (\beta_{x}^{(1)} = 1, \beta_{x}^{(0)} = 1)
% \end{equation} 
% The former substructure has been suggested by \cite{Haberman2011a} as a simpler structure to resolve some stability issues of the original model. The latter substructure is usually referred to as the Age-Period-Cohort (APC) model and has a long-standing tradition in the fields of medicine and demography (see, e,g,, \cite{Clayton1987} and \cite{Hobcraft1982}), but has not been widely used in the actuarial literature until considered by \cite{Currie2006}. Mortality projections for these models are typically derived from time series forecast of the estimated $\kappa_t$ and $\gamma_{t-x}$, generated using univariate ARIMA processes under the assumption of independence between the period and the cohort effects. The \rpkg{ilc} package \citep{Butt2014} implements the Renshaw and Haberman model and its substructures together with the LC model under a Poisson regression framework.

%   \subsection{The Cairns-Blake-Dowd model}\label{sec:CBDReview}
%   
% In addition to the Lee-Carter method several alternative mortality projection approaches have been proposed. The two factor model introduced by  \cite{Cairns2006a} is perhaps one of the most prominent ones. The Cairns-Blake-Dowd (CBD) model assumes the following specification for the mortality rates $q_{xt}$:
% \begin{equation}\label{eq:CBD}
% \logit q_{xt} = \log\left(\frac{q_{xt}}{1-q_{xt}}\right) = \kappa_t^{(1)} + (x-\bar{x})\kappa_t^{(2)}
% \end{equation}
% where $\bar{x}$ is the average age in the data and $\kappa^{(1)}_t$ and $\kappa^{(2)}_t$ are period effects. Whereas the Lee-Carter method assumes a non-parametric structure for the age effects, the CBD model treats age as a continuous covariate which enters the model in a linear way on the logit scale. This characteristic makes the model relatively parsimonious but restricts its applicability to older ages: a linear function cannot capture younger ages features such as the accident hump or infant mortality. As opposed to the Lee-Carter model (which considers only one period factor), the CBD model includes two period factors. As a result, the model can capture the imperfect correlation structure in mortality rates improvements across ages.\\
% 
% \cite{Cairns2009} consider three extensions to the original CBD model by incorporating a quadratic age term and a cohort effect:
% \begin{align}
% \logit q_{xt} &= \kappa_t^{(1)} + (x-\bar{x})\kappa_t^{(2)} + \gamma_{t-x} \label{eq:M6}\\
% \logit q_{xt} &= \kappa_t^{(1)} + (x-\bar{x})\kappa_t^{(2)} + \kappa^{(3)}_t \left((x-\bar{x})^2-\hat{\sigma}^2_x\right) + \gamma_{t-x} \label{eq:M7}\\
% \logit q_{xt} &= \kappa_t^{(1)} + (x-\bar{x})\kappa_t^{(2)} + (x_c - x)\gamma_{t-x} \label{eq:M8}
% \end{align}
% where $\hat{\sigma}^2_x$ is the average value of $(x-\bar{x})^2$ and $x_c$ is a constant parameter to be estimated. The models defined in equations \eqref{eq:M6}, \eqref{eq:M7}, and \eqref{eq:M8} are typically referred to as models M6, M7, and M8, respectively. These three models along with the original CBD model, the Lee-Carter model (using Poisson maximum likelihood) and the APC model are implemented in the \rpkg{LifeMetrics} \R software. This package, which is not on CRAN, is available at  \url{http://www.macs.hw.ac.uk/~andrewc/lifemetrics/}.\\
% 
% \cite{Plat2009a} has combined the CBD method with some features of the Lee-Carter method to produce a model that is suitable for full age ranges, captures the cohort effect, and has an imperfect correlation structure. The proposed model is
% \begin{equation}\label{eq:Plat}
% \log m_{xt} = \alpha_x + \kappa^{(1)}_t + (\bar{x}-x)\kappa^{(2)}_t + (\bar{x}-x)^+\kappa^{(3)}_t + \gamma_{t-x}
% \end{equation}
% where $(\bar{x}-x)^+=\max(0, \bar{x}-x)$. The  $\alpha_x$ term is the same as in the Lee-Carter method, allowing the model to capture the shape of the mortality rates for the full age range. In the cases where only older ages are of interests, \cite{Plat2009a} suggests to drop the third period term from \eqref{eq:Plat} to have the model   
% \begin{equation}\label{eq:PlatOld}
% \log m_{xt} = \alpha_x + \kappa^{(1)}_t + (\bar{x}-x)\kappa^{(2)}_t + \gamma_{t-x}.
% \end{equation}
% This reduced Plat model is essentially the M6 model of \cite{Cairns2009} with an added static age term $\alpha_x$.\\
% 
% In order to provide a better fit to ages 5 to 20, \cite{OHare2012} have modified the Plat model by incorporating a quadratic lower age effect:
% \begin{equation}\label{eq:Ohare}
% \log m_{xt} = \alpha_x + \kappa^{(1)}_t + (\bar{x}-x)\kappa^{(2)}_t +  \left((\bar{x}-x)^+-\left[(\bar{x}-x)^+\right]^2\right)\kappa^{(3)}_t + \gamma_{t-x}
% \end{equation}

  

 \section{GAPC stochastic mortality models with \rpkg{StMoMo}}\label{sec:GAPC_StMOMO}
 
 \revI{The \rpkg{StMoMo} package provides an \R implementation of the GAPC family of stochastic mortality models using the standard \rcode{S3} object-oriented system. \rpkg{StMoMo} can be installed with the code:}
 <<installStMoMo, eval = FALSE>>=
install.packages("StMoMo")
@
\noindent 
 \revI{The package is loaded within \R as follows:}
<<loadStMoMo, eval=TRUE,echo = TRUE, message=FALSE, warning=FALSE>>=
library(StMoMo)
@

 
 \revI{In package \rpkg{StMoMo}, GAPC stochastic mortality models are constructed using the \rcode{StMoMo} function}. This function takes as input information on the \textit{link function} (and the associated distributional assumption), the \textit{predictor structure}, and the \textit{set of parameter constraints} to create an object of \vlad{the} type \rcode{"StMoMo"} representing a GAPC mortality model:
\begin{itemize}
 \item The argument \rcode{link} defines the \textit{link function} and \revI{\textit{the random component}} associated with the mortality model. Setting \rcode{link = "log"} assumes that deaths follow a Poisson distribution and uses a log link targeting the force of mortality $\mu_{xt}$,  while setting \rcode{link = "logit"} assumes that deaths follow a Binomial distribution and uses a logit link targeting one-year mortality rates $q_{xt}$.
 
 \item The \textit{predictor} of the model is defined via arguments \rcode{staticAgeFun}, \rcode{periodAgeFun} and \rcode{cohortAgeFun}.  Argument \rcode{staticAgeFun} is a logical \vlad{variable} indicating whether the model has a static age function $\alpha_x$ or not. Argument \rcode{periodAgeFun} is a list of length $N$ containing the definitions of the period age-modulating parameters $\beta_x^{(i)}$, $i=1,\ldots,N$, with each entry being \vlad{either} \rcode{"NP"} for non-parametric age terms, \rcode{"1"} for $\beta_x^{(i)}=1$, or a predefined parametric function of age\footnote{\revI{Note that we can define a model with no age-period terms ($N=0$) by making \rcode{periodAgeFun = NULL}.}}.  Argument \rcode{cohortAgeFun} defines the cohort age modulating parameter $\beta_x^{(0)}$ and can take values  \rcode{"NP"} for non-parametric age terms, \rcode{"1"} for $\beta_x^{(0)}=1$, a predefined parametric function of age, or \rcode{NULL} if the model does not have a cohort effect.
 
 \item The \textit{set of parameter constraints} are defined via \vlad{the} argument \rcode{constFun} which is a user-defined \revI{implementation of the \textit{constraint function} $v$}  mapping an arbitrary \revI{vector} of \vlad{parameters} to a \revI{vector} of \vlad{transformed parameters} satisfying the model constraints.
\end{itemize}
\revI{We note that due to limitations of the \R functions used for fitting \rcode{"StMoMo"} objects to
data (see Section \ref{sec:Fitting}), the current version \rpkg{StMoMo} does not support models combining parametric and non-parametric age-modulating functions, $\beta_x^{(i)}$, $0=1,\ldots,N$. However, such model are not typically considered and the majority of models proposed in the literature are either extension of the Lee-Carter model with all age-modulating terms being non-parametric or extension of the CBD model with all age-modulating terms being parametric.}\footnote{\revI{For instance, a model with predictor structure $\eta_{xt} = \alpha_x + (x-\bar{x})\kappa^{(1)}_t + \beta^{(2)}_x\kappa^{(1)}_t$, corresponding to \rcode{StMoMo(staticAgeFun = TRUE, periodAgeFun = c(f1, "NP"))}, with \rcode{f1 <- function(x, ages) x - mean(ages)}, is not supported.}}  \\

\renewcommand{\arraystretch}{1.5}
\begin{table}[!tb]
  \centering
    \caption{Model structures considered in this paper.}\label{tab:Models}		
		\begin{tabular}{m{2cm} m{10cm} }
		\toprule		
		 \textbf{Model}& \textbf{Predictor} \tabularnewline \midrule		
		LC           & $\eta_{xt} = \alpha_x + \beta^{(1)}_x \kappa^{(1)}_t$ \tabularnewline 
		CBD                   & $\eta_{xt} = \kappa_t^{(1)} + (x-\bar{x})\kappa_t^{(2)}$ \tabularnewline  	
    APC                  & $\eta_{xt} = \alpha_x + \kappa^{(1)}_t+ \gamma_{t-x}$ \tabularnewline
		RH & $\eta_{xt} = \alpha_x + \beta^{(1)}_x \kappa^{(1)}_t + \gamma_{t-x}$\tabularnewline  	
		M7                   & $\eta_{xt} = \kappa_t^{(1)} + (x-\bar{x})\kappa_t^{(2)} + \left((x-\bar{x})^2-\hat{\sigma}_x^2\right)\kappa_t^{(3)} + \gamma_{t-x}$ \tabularnewline
    PLAT                   & $\eta_{xt} = \alpha_x + \kappa_t^{(1)} + (\bar{x}-x)\kappa_t^{(2)} + \gamma_{t-x}$ \tabularnewline
		\bottomrule
		\end{tabular}		
\end{table}
\renewcommand{\arraystretch}{1}

In order to illustrate the creation of particular GAPC mortality models and other \vlad{capabilities} of \rpkg{StMoMo}, in the rest of this paper we will focus on the models summarised in Table \ref{tab:Models}. From now onwards, LC stands for the Lee-Carter model; CBD for the original Cairns-Blake-Dowd model; APC for the Age-Period-Cohort model of \citet{Currie2006}; RH for the cohort extension of the Lee-Carter model defined in equation \eqref{eq:CohortEffectH1} and proposed by \citet{Renshaw2006}; M7 for the quadratic CBD model defined in equation \eqref{eq:M7}; and PLAT for the reduced Plat model defined previously in equation \eqref{eq:PlatOld}. For the sake of comparability, in all cases we will assume a Binomial distribution of deaths and use the logit function to link $q_{xt}$ to the predictor structure $\eta_{xt}$. \\

Below, we show how to define \vlad{each of the models in Table \ref{tab:Models} using the package \rpkg{StMoMo}}.\\

\noindent
\textbf{Lee-Carter model.} The LC model under a Binomial setting can be defined using the following code:

<<defineLClong, eval = FALSE>>=
constLC <- function(ax, bx, kt, b0x, gc, wxt, ages){
  c1 <- mean(kt[1, ], na.rm = TRUE)
  c2 <- sum(bx[, 1], na.rm = TRUE)
  list(ax = ax + c1 * bx, bx = bx / c2, kt =  c2 * (kt - c1))   
} 
LC <- StMoMo(link = "logit", staticAgeFun = TRUE, periodAgeFun = "NP", 
             constFun = constLC)
@
\noindent
Recalling Section \ref{ex:LC}, we note that the constraint function \rcode{constLC} is the \rcode{R} implementation of transformation \eqref{eq:LCTransformations} with constants $c_1$ and $c_2$ calculated using Equation \eqref{eq:LCConstants} to impose the constraints defined in equation \eqref{eq:LCConstraints}. The \rpkg{StMoMo} package also contains the function \rcode{lc} to facilitate the definition of Lee-Carter models. Hence, we could define the LC model using the much simpler predefined command:  
<<defineLC, eval=TRUE, cache=TRUE>>=
LC <- lc(link = "logit")
@

\noindent
\textbf{CBD model.} To define the CBD model we use the following commands:
<<defineCBDlong, eval = FALSE>>=
f2 <- function(x, ages) x - mean(ages)
CBD <- StMoMo(link = "logit", staticAgeFun = FALSE, 
              periodAgeFun = c("1", f2))
@
\noindent
Here, we note that function \rcode{f2} defines the second age-modulating parameter $\beta^{(2)}_x = x - \bar{x}$ and that a \rcode{constFUN} argument needs not be provided since the CBD model does not have identifiability issues. Alternatively, we can define the CBD model using the predefined function \rcode{cbd}:
<<defineCBD, eval=TRUE, cache=TRUE>>=
CBD <- cbd()
@

\noindent
\textbf{APC model, RH model and model M7.} These three models could be defined by implementing explicitly the discussions in Sections \ref{ex:RH}, \ref{ex:APC} and \ref{ex:M7}. However, \rpkg{StMoMo} includes predefined functions \rcode{apc}, \rcode{rh}, \rcode{m7} that facilitate the definition of the APC model, the RH model and model M7, respectively\footnote{The \rpkg{StMoMo} package also includes functions \rcode{m6} and \rcode{m8} for defining models M6 and M8. We also note that the \citet{Renshaw2006} cohort extension of the Lee-Carter in Equation \eqref{eq:LCCohortEffect} can be defined using function \rcode{rh} with argument \rcode{cohortAgeFun = "NP"}.}. Thus, these models are defined with the code:
<<defineAPC_RH_M7, eval=TRUE, cache=TRUE>>=
RH <- rh(link = "logit",  cohortAgeFun = "1")
APC <- apc(link = "logit")
M7 <- m7()
@  

\noindent
\textbf{PLAT model.} Package \rpkg{StMoMo} does not include a predefined function for the Plat model. Nevertheless, recalling Section \ref{ex:Plat}, we can define the reduced Plat model using the code:
<<definePlat, eval=TRUE, cache=TRUE>>=
f2 <- function(x, ages) mean(ages) - x
constPlat <- function(ax, bx, kt, b0x, gc, wxt, ages){
  nYears <- dim(wxt)[2]
  x <- ages
  t <- 1:nYears
  c <- (1 - tail(ages, 1)):(nYears - ages[1])
  xbar <- mean(x)
  #\sum g(c)=0, \sum cg(c)=0, \sum c^2g(c)=0
  phiReg <- lm(gc ~ 1 + c + I(c^2), na.action = na.omit)
  phi <- coef(phiReg)
  gc <- gc - phi[1] - phi[2] * c - phi[3] * c^2 
  kt[2, ] <- kt[2, ] + 2 * phi[3] * t
  kt[1, ] <- kt[1, ] + phi[2] * t + phi[3] * (t^2 - 2 * xbar * t)  
  ax <- ax + phi[1] - phi[2] * x + phi[3] * x^2
  #\sum kt[i, ] = 0
  ci <- rowMeans(kt, na.rm = TRUE)
  ax <- ax + ci[1] + ci[2] * (xbar - x) 
  kt[1, ] <- kt[1, ] - ci[1]
  kt[2, ] <- kt[2, ] - ci[2]
  list(ax = ax, bx = bx, kt = kt, b0x = b0x, gc = gc)
}
PLAT <- StMoMo(link = "logit", staticAgeFun = TRUE,
                  periodAgeFun = c("1", f2), cohortAgeFun = "1",
                  constFun = constPlat)
@
\noindent
We note that the constraint function \rcode{constPlat} is the \rcode{R} implementation of transformations \eqref{eq:PlatTrans2} and \eqref{eq:PlatTrans1} omitting the terms involving $\kappa^{(3)}$ and $c_3$, and with constants $\phi_1$, $\phi_2$, $\phi_3$ obtained via the linear regression defined in \eqref{eq:PlatConstFun1} and constant $c_1$ and $c_2$ as in Equation \eqref{eq:PlatConstFun3}.  Function \rcode{constPlat} imposes the constraints in Equation \eqref{eq:PlatConst}.\\

\section{Model fitting}\label{sec:Fitting}
Parameter estimates of GAPC stochastic mortality models can be obtained by maximising the model log-likelihood, which is given by
\begin{equation*}\label{eq:PoissonLogL}
\LogL({d_{xt}},{\hat{d}_{xt}}) = \sum_{x}  \sum_{t}\omega_{xt}\left\{{d_{xt}}\log{\hat{d}_{xt}}-\hat{d}_{xt}-\log{d_{xt}!}\right\}
\end{equation*}
in the case of a Poisson distribution of deaths, and by
\begin{equation*}\label{eq:BinomialLogL}
\LogL({d_{xt}},{\hat{d}_{xt}}) = \sum_{x}  \sum_{t}\omega_{xt}\left\{{d_{xt}}\log\left(\frac{\hat{d}_{xt}}{E^0_{xt}}\right) + (E^0_{xt}-d_{xt})\log\left(\frac{E^0_{xt}-\hat{d}_{xt}}{E^0_{xt}}\right) + \binom{E^0_{xt}}{d_{xt}} \right\}
\end{equation*}
in the case of a Binomial distribution of deaths. \vlad{In both cases, $\omega_{xt}$ are weights taking the value $0$ if a particular $(x,t)$ data cell is omitted or $1$ if the cell is included, and}
\[\hat{d}_{xt}=E_{xt}\,g^{-1}\left(\alpha_x + \sum_{i=1}^N \beta_x^{(i)}\kappa_t^{(i)} + \beta_x^{(0)}\gamma_{t-x}\right)\] is the expected number of deaths predicted by the model, with $g^{-1}$ denoting the inverse of the \textit{link function} $g$.\\

In the mortality literature, maximisation of the log-likelihood is typically performed using the Newton-Raphson iterative procedure tailored for each model (see e.g. \citet{Brouhns2002}, \citet{Renshaw2006} and  \citet{Cairns2009}). This is in fact the approach implemented in \vlad{the} packages \rpkg{ilc} and \rpkg{LifeMetrics}. Nonetheless, as discussed extensively by \cite{Currie2014}, many stochastic mortality models are examples of generalised linear models or generalised non-linear model, which facilitates their fitting using standard statistical software\footnote{\citet{Haberman2011a} have also noticed this fact and profit from GLM facilities in standard statistical packages when fitting CBD type models.}. \citet{Currie2014} exemplifies this fact by fitting several stochastic mortality models in \R using the standard function \rcode{glm} or the function \rcode{gnm} of \vlad{the} package \rpkg{gnm} \citep{Turner2012}\footnote{\revI{\citet{Debon2010} also discuss the use of the package \rpkg{gnm} for fitting Lee-Carter type models.}}.\\

\rpkg{StMoMo} provides the generic function \rcode{fit} for estimating the parameters of GAPC mortality models. In line with the remarks of \citet{Currie2014}, the corresponding \rcode{S3} method for objects of \vlad{the} class \rcode{"StMoMo"} relies on function \rcode{gnm} of \vlad{the} package \rpkg{gnm} to estimate the parameters of a GAPC model\footnote{The current version of \rpkg{StMoMo} only includes the \rcode{S3} method \rcode{fit.StMoMo} for objects of class \rcode{"StMoMo"}, but future versions may provide tailored \rcode{fit} methods for particular stochastic mortality models.}. Internally, this is accomplished by constructing the equivalent \rpkg{gnm} formulation of the GAPC mortality model\footnote{We note that when all the $\beta^{(i)}_x$ are parametric functions of age, the model is a GLM and therefore \rcode{gnm} by default resorts to the \rcode{glm} function of \rcode{R} when fitting the parameter of the model.}
.  For instance, the \rpkg{gnm} formula of the Binomial \rcode{LC} model created before is
<<showGNMformulaLC, eval=TRUE, cache=TRUE, dependson="defineLC", linewidth = 70>>=
LC$gnmFormula
@  
\noindent
while the \rpkg{gnm} formula of the Binomial \rcode{CBD} model defined before is
<<showGNMformulaCBD, eval=TRUE, cache=TRUE, dependson="defineCBD", linewidth = 70>>=
CBD$gnmFormula
@  

% \begin{figure}[tb]  
% \centering
% <<plotEWdata, eval=TRUE,  fig.width=5, fig.height=4,out.width="7.5cm", results='hide', warning=FALSE, message=FALSE, cache=TRUE, echo=FALSE>>=
% matplot(EWMaleData$ages, log(EWMaleData$Dxt / EWMaleData$Ext), type = "l", 
%         lty = 1, col = rainbow(60), xlab = "age", ylab = "log death rates", 
%         main = "EW: male death rates (1961-2011)")
% @
% \caption{England and Wales male mortality data for 1961-2011. In the graph warmer colours indicate earlier years while colder colours represent more recent years.}\label{fig:EWData}  
% \end{figure}

We now illustrate the usage of the function \rcode{fit} of \vlad{the} package \rpkg{StMoMo} by fitting the six models defined before to England and Wales mortality data. The object \rcode{EWMaleData} included in \vlad{the} package \rpkg{StMoMo} contains data on deaths and central exposures for England and Wales males for the period 1961-2011 and for ages 0-100 obtained from the \cite{HumanMortalityDatabase2014}. %The code below produces a graphical summary (see Figure \ref{fig:EWData}) of these data\footnote{Data for other countries from the \cite{HumanMortalityDatabase2014} can be easily downloaded in the appropriate format using the function \rcode{hmd.mx} of package \rpkg{demography}.}: 
% <<plotEWdatashow, eval=FALSE>>=
% matplot(EWMaleData$ages, log(EWMaleData$Dxt / EWMaleData$Ext), type = "l", 
%         lty = 1, col = rainbow(60), xlab = "age", ylab = "log death rates", 
%         main = "EW: male death rates (1961-2011)")
% @
% \noindent
\pietro{However, in our examples we concentrate} on ages 55 to 89 as the CBD model and the M7 model have been particularly designed to fit higher ages. Additionally, since some models include cohort effects and in agreement with the usual practice (see e.g \citet{Cairns2009} and \citet{Haberman2011a}), we exclude (by setting $\omega_{xt}=0$) all cohorts that have fewer than three observations.\\

Models LC, APC, CBD, M7 and PLAT can be fitted to England and Wales male mortality data for ages 55 to 89 using the code:
<<fit1, eval=TRUE, cache=TRUE, results='hide', dependson=c("defineLC", "defineAPC_RH_M7", "definePlat") >>=
Dxt <- EWMaleData$Dxt
Ext <- EWMaleData$Ext + 0.5 * EWMaleData$Dxt
ages <- EWMaleData$ages
years <- EWMaleData$years
ages.fit <- 55:89
wxt <- genWeightMat(ages = ages.fit, years = years, clip = 3)
LCfit <- fit(LC, Dxt = Dxt, Ext = Ext, ages = ages, years = years, 
             ages.fit = ages.fit, wxt = wxt)
APCfit <- fit(APC, Dxt = Dxt, Ext = Ext, ages = ages, years = years, 
             ages.fit = ages.fit, wxt = wxt)
CBDfit <- fit(CBD, Dxt = Dxt, Ext = Ext, ages = ages, years = years, 
             ages.fit = ages.fit, wxt = wxt)
M7fit <- fit(M7, Dxt = Dxt, Ext = Ext, ages = ages, years = years, 
             ages.fit = ages.fit, wxt = wxt)
PLATfit <- fit(PLAT, Dxt = Dxt, Ext = Ext, ages = ages, years = years, 
             ages.fit = ages.fit, wxt = wxt)
@
\noindent
From this code we note the following:
\begin{itemize}
  \item In order to match the logit-Binomial setting used before in the definition of the mortality models, initial exposures are approximated by transforming the available central exposures.
  \item The first and last three cohorts years are excluded from the fitting via argument \rcode{wxt}. The appropriate 0-1 weighting matrix, \rcode{wxt}, is constructed using the utility function \rcode{genWeightMat} of package \rpkg{StMoMo}.
\end{itemize}

The fitting of the RH model requires some care as it is well known that fitting cohort extensions of the Lee-Carter models is problematic \citep{Hunt2014a}. In particular, \cite{Currie2014} has encountered convergence issues when using package \rpkg{gnm} to fit the RH model. As a possible way to circumvent these issues, \cite{Currie2014} suggests the use of appropriate starting values when fitting model RH. This can be achieved in function \rcode{fit} via input arguments \rcode{start.ax}, \rcode{start.bx}, \rcode{start.kt}, \rcode{start.b0x}, and \rcode{start.gc}. Using the parameters of the Lee-Carter model as starting values, model RH can be fitted with the code:
  
<<fitRH, eval=TRUE, cache=TRUE, results='hide', dependson="defineAPC_RH_M7">>=
RHfit <- fit(RH, Dxt = Dxt, Ext = Ext, ages = ages, years = years, 
             ages.fit = ages.fit, wxt = wxt, start.ax = LCfit$ax, 
             start.bx = LCfit$bx, start.kt = LCfit$kt)
@

The output from the function \rcode{fit} is an object of \vlad{the} class \rcode{"fitStMoMo"} including, among other things, the following information:
\begin{itemize}  
  \item \rcode{model}: the \rcode{"StMoMo"} object defining the underlying GAPC stochastic mortality model;
  \item \rcode{ax}, \rcode{bx}, \rcode{kt}, \rcode{b0x}, \rcode{gc}: the estimated parameters;
  \item \rcode{loglik}: the log-likelihood of the model;
  \item \rcode{deviance}: the model deviance;
  \item \rcode{nobs}: the number of observations in the data;
  \item \rcode{npar}: the effective number of parameters of the model;  
  \item \rcode{fittingModel}: the output of the \rcode{gnm} call used to fit the model.
\end{itemize}
There are {\rcode{print}}, \rcode{plot}, \rcode{fitted}, and \rcode{residuals} methods for the \rcode{"fitStMoMo"} class. For instance, Figures \ref{fig:ParamLC}, \ref{fig:ParamCBD} and \ref{fig:ParamAPC} depicting the fitted parameters of the LC model, the CBD model and the APC model, respectively, were produced with the code:
<<plotParam, eval=FALSE>>=
plot(LCfit, nCol = 3)
plot(CBDfit, parametricbx = FALSE)
plot(APCfit, parametricbx = FALSE, nCol = 3)
@ 

\begin{figure}[!bt]  
\centering
<<plotLC, fig.width=8, fig.height = 2.8, out.width="16cm", echo=FALSE, cache=TRUE, dependson="fit1">>=
plot(LCfit, nCol = 3)
@ 
\caption{Parameters for the Lee-Carter (LC) model fitted to the England and Wales male population for ages 55-89 and the period 1961-2011.}\label{fig:ParamLC}  
\end{figure}

\begin{figure}[!tb]  
\centering
<<plotCBD, fig.width=8, fig.height = 4.2, out.width="11cm", echo=FALSE, cache=TRUE, dependson="fit1">>=
plot(CBDfit, parametricbx = FALSE)
@ 
\caption{Parameters for the CBD model fitted to the England and Wales male population for ages 55-89 and the period 1961-2011.}\label{fig:ParamCBD}  
\end{figure}

\begin{figure}[!bt]  
\centering
<<plotAPC, fig.width=8, fig.height = 2.8, out.width="16cm", echo=FALSE, cache=TRUE, dependson="fit1">>=
plot(APCfit, parametricbx = FALSE, nCol = 3)
@ 
\caption{Parameters for the APC model fitted to the England and Wales male population for ages 55-89 and the period 1961-2011.}\label{fig:ParamAPC}  
\end{figure}

%\begin{itemize}  
% 	 \draft{When fitting M7* \cite{Haberman2011a} employ a two-stage fitting strategy estimating first $\alpha_x$ as $\hat{\alpha}_x = \sum_t \log\left(\slfrac{d_{xt}}{e_{xt}}\right)$ and the treating $\hat{\alpha}_x$ as an offset when estimating the other parameters . This can be achieved using parameter \rcode{oxt} in \rcode{fit} using the code:}

% <<fitM7star, eval=TRUE, cache=TRUE, results='hide', dependson="defineM7star">>=
% axhat <-rowMeans(log(Dxt / Ext))[which (ages %in% ages.fit)]
% M7starfit <- fit(M7star, Dxt = Dxt, Ext = Ext, ages = ages, years = years, 
%              ages.fit = ages.fit, wxt = wxt, oxt = axhat)
% @ 
% 
% \item \draft{Example of summary}
%  
% 
% \end{itemize}


\section{Goodness-of-fit analysis}\label{sec:GoF}


<<resAll, eval=TRUE, echo=FALSE, cache=TRUE, dependson=c("fit1","fitRH")>>=
LCres <- residuals(LCfit)
CBDres <- residuals(CBDfit)
APCres <- residuals(APCfit)
RHres <- residuals(RHfit)
M7res <- residuals(M7fit)
PLATres <- residuals(PLATfit)
@ 
\pietro{The goodness-of-fit of mortality models is typically analysed} by inspecting the residuals of the fitted model. Regular patterns in the residuals indicate the inability of the model to describe all the features of the data appropriately. With a Poisson or Binomial random component, it is appropriate to look at the scaled deviance residuals defined as:  

\begin{equation*}\label{eq:Residuals}
r_{xt}=\mathrm{sign}(d_{xt}-\hat{d}_{xt})\sqrt{\frac{\mathrm{dev}(x,t)}{\hat{\phi}}}, \qquad \hat{\phi} = \frac{D(d_{xt},\hat{d}_{xt})}{K-\nu},
\end{equation*}
where 
\[\mathrm{dev}(x,t) = 2d_{xt}\log\left(\frac{d_{xt}}{\hat{d}_{xt}}\right) - (d_{xt} - \hat{d}_{xt})\]
for a Poisson random component, or
\[\mathrm{dev}(x,t) = 2d_{xt}\log\left(\frac{d_{xt}}{\hat{d}_{xt}}\right) + (E^0_{xt}-d_{xt})\log\left(\frac{E^0_{xt}-d_{xt}}{E^0_{xt}-\hat{d}_{xt}}\right)\]
for a Binomial random component. \pietro{Further,
\[D(d_{xt},\hat{d}_{xt}) = \sum_{x}  \sum_{t}\omega_{xt}\mathrm{dev}(x,t)\] 
is the total deviance of the model, $K=\sum_x\sum_t\omega_{xt}$ is the number of observations in the data and} $\nu$ is the effective number of parameters in the model.\\

In \rpkg{StMoMo} standardised deviance residuals can be obtained with the generic function \rcode{residuals} applied to a fitted stochastic mortality model of \vlad{the} class \rcode{"fitStMomo"}. For example, to obtain the residuals of the LC model and of the CBD model fitted before we use the commands:
<<resLC_CBD, eval=FALSE>>=
LCres <- residuals(LCfit)
CBDres <- residuals(CBDfit)
@ 
\noindent
Graphs of these residuals can be produced using the generic function \rcode{plot}. This function supports, via argument \rcode{type}, three types of plots: 
\begin{itemize}
\item scatter plots of residuals by age, period and cohort such as those extensively used in \citet{Haberman2011a};
\item black and white sign-plots of the residuals such as those used in \citet{Cairns2009} and \citet{Lovasz2011}; and 
\item colour maps (heat-maps) of the residuals. 
\end{itemize}


Figure \ref{fig:ResHeatmap} presents heat-maps of the deviance residuals for the six models fitted to the England and Wales male mortality experience. These charts were produced using function \rcode{plot} with option \rcode{type = "colourmap"}. For instance, Figure \ref{fig:ResHeatmapCBD} was obtained with the code:  
<<plotCBDresShow, eval=FALSE>>=
plot(CBDres, type = "colourmap", reslim = c(-3.5, 3.5))
@ 
\noindent
From Figure \ref{fig:ResHeatmap} we see that models LC, CBD and APC display strong residual patterns while the residuals of models RH, M7 and PLAT look reasonably random. The APC model shows a strong clustering of residuals due to its inability to allow for varying improvement rates with age. The LC and CBD models, which do not incorporate a cohort effect, show very marked diagonals patterns indicating the inability of these models to capture the well-known cohort effect observed in the England and Wales population \citep{Willets2004}. The issues with the fit of the LC and CBD models become more evident when looking at scatter plots of the residuals by age, period and cohort. Such plots for the LC model (Figure \ref{fig:ResScatterLC}) and the CBD model (Figure \ref{fig:ResScatterCBD}) can be produced using argument \rcode{type = "scatter"} of function \rcode{plot} via the commands:
<<plotResShow, eval=FALSE>>=
plot(LCres, type = "scatter", reslim = c(-3.5, 3.5))
plot(CBDres, type = "scatter", reslim = c(-3.5, 3.5))
@ 
\noindent
The right panels in Figure \ref{fig:ResScatter} clearly show that the LC and CBD models are unable to capture the cohort effect. In addition, the left panel in Figure \ref{fig:ResScatterCBD} reveals some strong patterns by age, reflecting the the lack of a quadratic age term in the CBD which may be necessary to capture the commonly observed curvature of the mortality rates in a logit scale.\\


\begin{figure}[!htb]  
    \centering
    \subfloat[LC]{
    <<plotLCres, fig.width=4, fig.height=3,out.width="7.8cm", echo=FALSE, cache=TRUE, dependson="resAll">>=
par(mar=c(4.5,4,1,1))
plot(LCres, type = "colourmap", reslim = c(-3.5,3.5))
@ 
}
  	\subfloat[CBD]{
    <<plotCBDres, fig.width=4, fig.height=3,out.width="7.8cm", echo=FALSE, cache=TRUE, dependson="resAll">>=
par(mar=c(4.5,4,1,1))
plot(CBDres, type = "colourmap", reslim = c(-3.5,3.5))
@ 
\label{fig:ResHeatmapCBD}}\,		

\subfloat[APC]{
    <<plotAPCres, fig.width=4, fig.height=3,out.width="7.8cm", echo=FALSE, cache=TRUE, dependson="resAll">>=
par(mar=c(4.5,4,1,1))
plot(APCres, type = "colourmap", reslim = c(-3.5,3.5))
@ 
}
\subfloat[RH]{
    <<plotRHres, fig.width=4, fig.height=3,out.width="7.8cm", echo=FALSE, cache=TRUE, dependson="resAll">>=
par(mar=c(4.5,4,1,1))
plot(RHres, type = "colourmap", reslim = c(-3.5,3.5))
@ 
}\,
\subfloat[M7]{
    <<plotM7res, fig.width=4, fig.height=3,out.width="7.8cm", echo=FALSE, cache=TRUE, dependson="resAll">>=
par(mar=c(4.5,4,1,1))
plot(M7res, type = "colourmap", reslim = c(-3.5,3.5))
@ 
}
\subfloat[PLAT]{
    <<plotPLATres, fig.width=4, fig.height=3,out.width="7.8cm", echo=FALSE, cache=TRUE, dependson="resAll">>=
par(mar=c(4.5,4,1,1))
plot(PLATres, type = "colourmap", reslim = c(-3.5,3.5))
@ 
}
		\caption{Heat-maps of deviance residuals for different model fitted to the England and Wales males population for ages 55-89 and the period 1961-2011.}\label{fig:ResHeatmap}		
\end{figure}


\begin{figure}[!htb]  
    \centering
    
    \subfloat[LC]{
    <<plotLCresScatter, fig.width=10, fig.height=3,out.width="16cm", echo=FALSE, cache=TRUE, dependson="resAll">>=
plot(LCres, type = "scatter", reslim = c(-3.5,3.5))
@ 
\label{fig:ResScatterLC}}\,
\subfloat[CBD]{
    <<plotCBDresScatter, fig.width=10, fig.height=3,out.width="16cm", echo=FALSE, cache=TRUE, dependson="resAll">>=
plot(CBDres, type = "scatter", reslim = c(-3.5,3.5))
@ 
\label{fig:ResScatterCBD}} 	
\caption{Scatter plots of deviance residuals for models LC and CBD fitted to the England and Wales males population for ages 55-89 and the period 1961-2011.}\label{fig:ResScatter}		
\end{figure}

When evaluating the goodness-of-fit of different models, it is generally anticipated that models with more parameters provide a better fit to the data. To rule out the possibility that the better fit observed in a model is the result of over-parametrisation and compare the relative performance of several models, it has become common in the mortality literature to use information criteria which modify the maximum likelihood criterion by penalising models with more parameters\footnote{For examples of such analysis see \citet[Section 6.1.1]{Cairns2009}, \citet[Section 3.3]{Haberman2011a}, and \citet[Section 4.1]{Lovasz2011}.}. Two of these criteria are the Akaike Information Criteria (AIC) and the Bayesian Information Criteria (BIC), defined as $AIC = 2\nu-2\LogL$ and $BIC = \nu \log K - 2\LogL$, respectively, with a lower value of AIC and BIC being preferable. In \R these information criteria can be computed using the generic functions \rcode{AIC} and \rcode{BIC}. For example, we can get the AIC and BIC of the CBD model as follows:
<<BICCBD, cache=TRUE, dependson="fit1">>=
AIC(CBDfit)
BIC(CBDfit)
@ 
<<BICtable, echo=FALSE, cache=TRUE, dependson="fit1", results='asis', message=FALSE, warning=FALSE>>=
library(xtable)
BICdf <-data.frame(list(LC = c(LCfit$npar[1], AIC(LCfit), BIC(LCfit)), 
                        CBD = c(CBDfit$npar[1], AIC(CBDfit), BIC(CBDfit)), 
                        APC = c(APCfit$npar[1], AIC(APCfit), BIC(APCfit)), 
                        RH = c(RHfit$npar[1], AIC(RHfit), BIC(RHfit)), 
                        M7 = c(M7fit$npar[1], AIC(M7fit), BIC(M7fit)), 
                        PLAT = c(PLATfit$npar[1], AIC(PLATfit), BIC(PLATfit))))
rownames(BICdf) <- c("Number of Parameters", "AIC", "BIC")
xtable(BICdf, dec = 1, caption = "Number of parameters, AIC and BIC values for different model fitted to the England and Wales males population for ages 55-89 and the period 1961-2011", center = "centering", file = "", floating = TRUE,label="tab:InfoCriteria",align="lcccccc", digits = 0)
@ 
\noindent
Table \ref{tab:InfoCriteria} presents AIC and BIC values for the six models fitted to the England and Wales male data. We note that both criteria lead to the same ranking of models with M7, PLAT, and RH being the best performing models. Overall, these results are consistent with the existing literature comparing single population models, where the Renshaw-Haberman extension of the Lee-Carter model and the M7 model have been identified as good candidates for modelling mortality in the England and Wales population \citep{Cairns2009,Haberman2011a}.   

\section{Forecasting and simulation with stochastic mortality models}\label{sec:Forecasting}
In the family of GAPC stochastic mortality models the dynamics of mortality are driven by the period indexes $\kappa^{(i)}_t$, $i=1,\ldots, N$, and the cohort index $\gamma_{t-x}$. Therefore, the forecasting and simulation of mortality rates requires the modelling of these indexes using time series techniques.\\

For the period indexes we follow the standard approach \citep{Cairns2006a, Cairns2011a, Haberman2011a, Lovasz2011} and use a multivariate random walk with drift. Specifically, we assume that

\begin{equation}\label{eq:GAPCPeriodIndex}
\boldsymbol{\kappa}_t = \boldsymbol{\delta} + \boldsymbol{\kappa}_{t-1} + \boldsymbol{\xi}_t^\kappa, \qquad \boldsymbol{\kappa}_t = 
\left(\!
    \begin{array}{c}
      \kappa_t^{(1)} \\
      \vdots\\
  		\kappa_t^{(N)}
    \end{array}
  \!\right), \qquad \boldsymbol{\xi}_t^\kappa\sim N(\mathbf{0},\Sigma),
\end{equation}
where \pietro{$\boldsymbol{\delta}$} is an $N$-dimension vector of drift parameters and $\Sigma$ is the $N\times N$ variance-covariance matrix of the multivariate white noise $\boldsymbol{\xi}_t^\kappa$.\\

As pointed out by \citet{Currie2014}, the main challenge when forecasting stochastic mortality models is specifying the dynamics of the cohort effect. To have a simple starting point, we follow previous studies \citep{Renshaw2006, Cairns2011a, Lovasz2011} and assume that the cohort index, $\gamma_{t-x}$, follows a univariate ARIMA process which is independent of the the period index, $\boldsymbol{\kappa}_t$. In general, we assume that $\gamma_c\equiv \gamma_{t-x}$ follows an ARIMA$(p,q,d)$ with drift, so that 
\begin{equation}\label{eq:GAPCCohortIndex}
\Delta^d\gamma_c =  \delta_0 + \phi_1\Delta^d\gamma_{c-1} +\cdots+ \phi_p\Delta^d\gamma_{c-p}+\epsilon_{c}+\delta_1\epsilon_{c-1}+\cdots+\delta_q\epsilon_{c-q},
\end{equation}
where $\Delta$ is the difference operator, $\delta_0$ is the drift parameter, $\phi_1,\dots,\phi_p$ are the autoregressive coefficients with $\phi_p\neq 0$,  $\delta_1,\dots,\delta_q$ are the moving average coefficients with $\delta_q\neq 0$ and $\epsilon_c$ is a Gaussian white noise process with variance $\sigma_{\epsilon}$.\\ 

The time series models in \eqref{eq:GAPCPeriodIndex} and \eqref{eq:GAPCCohortIndex} can be used to obtain projected (simulated) values of the period index  $\dot{\boldsymbol{\kappa}}_{t_n+s}:=\left(\dot{\kappa}^{(1)}_{t_n+s}, \ldots, \dot{\kappa}^{(N)}_{t_n+s}\right)'$ and cohort index $\dot{\gamma}_{t_n + s - x_1}$, $s=1,\ldots,h$, respectively, to derive forecasted (simulated)  values of the predictor
\begin{equation*}\label{eq:ProjectedPredictor}
\dot{\eta}_{x,t_n+s} = \alpha_x + \sum_{i=1}^N \beta_x^{(i)}\dot{\kappa}_{t_n+s}^{(i)} + \beta_x^{(0)}\dot{\gamma}_{t_n+s-x},
\end{equation*}
which can in turn be used to obtain forecasted (simulated) age-specific central mortality rates, $\dot{\mu}_{x,t_n+s}$ or age-specific one-year mortality rates, $\dot{q}_{x,t_n+s}$.\\

\renewcommand{\arraystretch}{1.5}
\begin{table}[tb]
  \centering
    \caption{ARIMA Models for the cohort effect for models APC, RH, M7 and PLAT.}\label{tab:ARIMAModels}  	
		\begin{tabular}{m{4cm} m{8cm} }
		\toprule		
		 \textbf{Mortality Model}& \textbf{Model for $\gamma_{t-x}$} \tabularnewline \midrule		
    APC                      & ARIMA$(1,1,0)$ with drift \tabularnewline
		RH                       & ARIMA$(1,1,0)$ with drift\tabularnewline  	
		M7                       & ARIMA$(2,0,0)$ with non-zero intercept\tabularnewline
    PLAT                     & ARIMA$(2,0,0)$ with non-zero intercept \tabularnewline
		\bottomrule
		\end{tabular}		
\end{table}
\renewcommand{\arraystretch}{1}


In \vlad{the} package \rpkg{StMoMo} the forecasting of GAPC stochastic mortality models is implemented via the generic method \rcode{forecast}. This function estimates and forecasts the multivariate random walk with drift in Equation \eqref{eq:GAPCPeriodIndex} using the approach described in \citet[Appendix B]{Haberman2011a}\footnote{This is implemented in function \rcode{mrwd} of \vlad{the} package \rpkg{StMoMo}.} and uses function \rcode{Arima} of package \rpkg{forecast} \citep{Hyndman2008} to estimate and forecast the ARIMA process of Equation \eqref{eq:GAPCCohortIndex}. For instance, if we assume that the cohort indexes of the APC, RH, M7, and PLAT model follow the ARIMA processes specified in Table \ref{tab:ARIMAModels}, 50-year ahead ($h=50$) central projections of the period indexes, cohort index, and one-year death probabilities for the England and Wales mortality experience can be obtain with the code:
    <<forAll, cache=TRUE, dependson=c("fit1","fitRH")>>=
LCfor <- forecast(LCfit, h = 50)
CBDfor <- forecast(CBDfit, h = 50)
APCfor <- forecast(APCfit, h = 50, gc.order = c(1, 1, 0))
RHfor <- forecast(RHfit, h = 50, gc.order = c(1, 1, 0))
M7for <- forecast(M7fit, h = 50, gc.order = c(2, 0, 0))
PLATfor <- forecast(PLATfit, h = 50, gc.order = c(2, 0, 0))
@ 
The output from the function \rcode{forecast} is an object of \vlad{the} class \rcode{"forStMoMo"} including, among other things, the following information:
\begin{itemize}  
  \item \rcode{rates}: a matrix with the central projection of the mortality rates, $\dot{\mu}_{x,t_n+s}$ or $\dot{q}_{x,t_n+s}$, $s=1,\ldots,h$;
  \item \rcode{kt.f}: a list containing information on the multivariate random walk with drift fitted to the period index $\boldsymbol{\kappa}_t$; and 
  \item \rcode{gc.f}: a list containing information on the ARIMA model fitted to the cohort index $\gamma_{t-x}$.
  \end{itemize}
There are {\rcode{print}} and \rcode{plot} methods for the \rcode{"forStMoMo"} class. For instance, plots of the forecast of the period indexes of models RH, M7 and PLAT (see Figure \ref{fig:PeriodIndexForecast}) can be produced using the code:  

<<plotForecastPeriodIndexShow, eval=FALSE>>=  
plot(RHfor, only.kt = TRUE)
plot(M7for, only.kt = TRUE, nCol = 3)
plot(PLATfor, only.kt = TRUE)
@

\noindent
Similarly, plots of the forecast of the cohort indexes of the APC, RH, M7 and PLAT models (see Figure \ref{fig:CohortIndexForecast}) can be obtained with the commands: 
<<plotForecastCohortIndexShow, eval=FALSE>>=  
plot(APCfor, only.gc = TRUE)
plot(RHfor, only.gc = TRUE)
plot(M7for, only.gc = TRUE)
plot(PLATfor, only.gc = TRUE)
@

Package \rpkg{StMoMo} also provides the function \rcode{simulate} for simulating trajectories from GAPC stochastic mortality models. To simulate the period index, $\boldsymbol{\kappa}_t$, \rpkg{StMoMo} implements a multivariate adaptation of Algorithm 2 in \citet{Haberman2009} without provision for parameter error\footnote{We note that Algorithm 2 in \citet{Haberman2009} is itself an adaptation of the prediction interval approach of \citet{Cairns2006a}.}, while to simulate the cohort index, $\gamma_{t-x}$, function \rcode{simulate} uses the equivalent \rcode{S3} method for objects of class \rcode{"Arima"} provided by \vlad{the} package \rpkg{forecast}. For example, the code below produces 500 simulated trajectories \revI{for the next 50 years} of the six stochastic mortality models fitted previously to the England and Wales male mortality experience:
<<simAll, cache=TRUE, dependson=c("fit1","fitRH")>>=
set.seed(1234)
nsim <- 500
LCsim <- simulate(LCfit, nsim = nsim, h = 50)
CBDsim <- simulate(CBDfit, nsim = nsim, h = 50)
APCsim <- simulate(APCfit, nsim = nsim, h = 50, gc.order= c(1, 1, 0))
RHsim <- simulate(RHfit, nsim = nsim, h = 50, gc.order= c(1, 1, 0))
M7sim <- simulate(M7fit, nsim = nsim, h = 50, gc.order= c(2, 0, 0))
PLATsim <- simulate(PLATfit, nsim = nsim, h = 50, gc.order= c(2, 0, 0))
@ 

\begin{landscape}
\begin{figure}[!htb]  
    \centering
    
    \subfloat[RH]{
    <<plotForecastPeriodIndexRH, fig.width=5.2, fig.height=4,out.width="7.5cm", echo=FALSE, cache=TRUE, dependson="forAll">>=
plot(RHfor, only.kt = TRUE, lwd = 2)
@ 
\label{fig:PeriodRH}}
\subfloat[PLAT]{
    <<plotForecastPeriodIndexPLAT, fig.width=10.4, fig.height=4,out.width="15cm", echo=FALSE, cache=TRUE, dependson="forAll">>=
plot(PLATfor, only.kt = TRUE, lwd = 2)
@ 
\label{fig:PeriodPLAT}}\,   
\subfloat[M7]{
    <<plotForecastPeriodIndexM7, fig.width=10.4, fig.height=2.8,out.width="22.5cm", echo=FALSE, cache=TRUE, dependson="forAll">>=
plot(M7for, only.kt = TRUE, nCol = 3, lwd = 2)
@ 
\label{fig:PeriodM7}}   
\caption{Forecast of the period indexes of the RH, M7 and PLAT models applied to the England and Wales males population for ages 55-89 and the period 1961-2011. Dashed lines represent central forecast and dotted lines represent 95\% prediction intervals.}\label{fig:PeriodIndexForecast}    
\end{figure}
\end{landscape}

\begin{figure}[!tb]  
    \centering    
    \subfloat[APC]{
    <<plotForecastCohortIndexAPC, fig.width=6, fig.height=4,out.width="8cm", echo=FALSE, cache=TRUE, dependson="forAll">>=
plot(APCfor, only.gc = TRUE, lwd = 2)
@ 
\label{fig:CohortAPC}}
    \subfloat[RH]{
    <<plotForecastCohortIndexRH, fig.width=6, fig.height=4,out.width="8cm", echo=FALSE, cache=TRUE, dependson="forAll">>=
plot(RHfor, only.gc = TRUE, lwd = 2)
@ 
\label{fig:CohortRH}}\,
    \subfloat[M7]{
    <<plotForecastCohortIndexM7, fig.width=6, fig.height=4,out.width="8cm", echo=FALSE, cache=TRUE, dependson="forAll">>=
plot(M7for, only.gc = TRUE, lwd = 2)
@ 
\label{fig:CohortM7}}
    \subfloat[PLAT]{
    <<plotForecastCohortIndexPLAT, fig.width=6, fig.height=4,out.width="8cm", echo=FALSE, cache=TRUE, dependson="forAll">>=
plot(PLATfor, only.gc = TRUE, lwd = 2)
@ 
\label{fig:CohortPLAT}}
\caption{Forecast of the cohort indexes of the APC, RH, M7 and PLAT models applied to the England and Wales males population for ages 55-89 and the period 1961-2011. Dashed lines represent central forecast and dotted lines represent 95\% prediction intervals.}\label{fig:CohortIndexForecast}    
\end{figure}
\noindent
The output from the function \rcode{simulate} is an object of \vlad{the} class \rcode{"simStMoMo"} including, among other things, the following information:
\begin{itemize}  
  \item \rcode{rates}: a three dimensional array with the future simulated mortality rates;
  \item \rcode{kt.s}: a list containing information on the simulated paths of the period index $\boldsymbol{\kappa}_t$; and 
  \item \rcode{gc.s}: a list containing information on the simulated paths of the cohort index $\gamma_{t-x}$.
\end{itemize}
This output can be used to extract sample trajectories from a model. For instance, Figure \ref{fig:SimRH}, which depicts 20 trajectories of the period index, cohort index and one-year death rates at age 65 from model RH, was produced with the code:  
<<plotSimulationRHShow, eval=FALSE>>=
#Plot period index
par(mfrow = c(1, 3))
plot(RHfit$years, RHfit$kt[1, ],
     xlim = range(RHfit$years, RHsim$kt.s$years),
     ylim = range(RHfit$kt, RHsim$kt.s$sim[1, , 1:20]), 
     type = "l", xlab = "year", ylab = "kt", main = "Period index")
matlines(RHsim$kt.s$years, RHsim$kt.s$sim[1, , 1:20], type = "l", lty = 1)
#Plot cohort index
plot(RHfit$cohorts, RHfit$gc, 
     xlim = range(RHfit$cohorts, RHsim$gc.s$cohorts),
     ylim = range(RHfit$gc, RHsim$gc.s$sim[, 1:20], na.rm = TRUE),
     type = "l", xlab = "year", ylab = "kt",
     main = "Cohort index (ARIMA(1,1,0) with drift)")
matlines(RHsim$gc.s$cohorts, RHsim$gc.s$sim[, 1:20], type = "l", lty = 1)
#Plot rates at age 65
qxt <- Dxt / Ext
plot(RHfit$years, qxt["65", ], xlim = range(RHfit$years, RHsim$years), 
     ylim = range(qxt["65", ], RHsim$rates["65", , 1:20]), type = "l", 
     xlab = "year", ylab = "rate", main = "Mortality rates at age 65")
matlines(RHsim$years, RHsim$rates["65", , 1:20], type = "l", lty = 1)
@
\begin{figure}[tb]  
\centering
<<plotSimulationRH, eval=TRUE,  fig.width=9, fig.height=3,out.width="16cm", results='hide', warning=FALSE, message=FALSE, cache=TRUE, echo=FALSE, dependson="simAll">>=
par(mfrow=c(1, 3))
plot(RHfit$years, RHfit$kt[1, ],
     xlim = range(RHfit$years, RHsim$kt.s$years),
     ylim = range(RHfit$kt, RHsim$kt.s$sim[1, , 1:20]), 
     type = "l", xlab = "year", ylab = "kt", main = "Period index")
matlines(RHsim$kt.s$years, RHsim$kt.s$sim[1, , 1:20], type = "l", lty = 1)
#Plot cohort index
plot(RHfit$cohorts, RHfit$gc, 
     xlim = range(RHfit$cohorts, RHsim$gc.s$cohorts),
     ylim = range(RHfit$gc, RHsim$gc.s$sim[, 1:20], na.rm = TRUE),
     type = "l", xlab = "year", ylab = "kt",
     main = "Cohort index (ARIMA(1,1,0) with drift)")
matlines(RHsim$gc.s$cohorts, RHsim$gc.s$sim[, 1:20], type = "l", lty = 1)
#Plot rates at age 65
qxt <- Dxt / Ext
plot(RHfit$years, qxt["65", ], xlim = range(RHfit$years, RHsim$years), 
     ylim = range(qxt["65", ], RHsim$rates["65", , 1:20]), type = "l", 
     xlab = "year", ylab = "rate", main = "Mortality rates at age 65")
matlines(RHsim$years, RHsim$rates["65", , 1:20], type = "l", lty = 1)
@
\caption{20 simulated trajectories of the period index $\kappa^{(1)}_t$, cohort index $\gamma_{t-x}$ and one-year death rates at age 65 $q_{xt}$ for model RH fitted to the England and Wales males population for ages 55-89 and the period 1961-2011.}\label{fig:SimRH}  
\end{figure}

We can also use the output from function \rcode{simulate} to produce fan charts depicting the uncertainty  associated with a model forecast. According to \citet{Cairns2011a}, such plots are central to the analysis of the plausibility of the forecast from a model, and can be used as a criterion when deciding upon what is the most appropriate model amongst a group of possible stochastic mortality models. Figure \ref{fig:Fancharts} shows fan charts depicting 50\%, 80\% and 95\% prediction intervals for mortality rates at ages 65, 75 and 85 for each of the six models fitted to the England and Wales experience. Figure \ref{fig:FanCBD}, for instance, was produced with the code:
<<plotCBDFanShow, eval=FALSE>>=
library(fanplot)
probs = c(2.5, 10, 25, 50, 75, 90, 97.5)
qxt <- Dxt / Ext 
#age 65
plot(CBDfit$years, qxt["65", ], xlim = c(1960, 2061), ylim = c(0.0025, 0.2),
     xlab = "year", ylab = "mortality rate (log scale)", 
     pch = 20, log = "y")
fan(t(CBDsim$rates["65", , ]), start = 2012, probs = probs, n.fan = 4,
    fan.col = colorRampPalette(c("black", "white")), ln = NULL)
#age 75
points(CBDfit$years, qxt["75", ], pch = 20)
fan(t(CBDsim$rates["75", , ]), start = 2012, probs = probs, n.fan = 4,
    fan.col = colorRampPalette(c("red", "white")), ln = NULL)
#age 75
points(CBDfit$years, qxt["85", ], pch = 20)
fan(t(CBDsim$rates["85", , ]), start = 2012, probs = probs, n.fan = 4,
    fan.col = colorRampPalette(c("blue", "white")), ln = NULL)
#labels
text(1965, qxt[c("65", "75", "85"), "1990"], 
     labels = c("x = 65", "x = 75", "x = 85"))
@ 
\noindent
From Figure \ref{fig:Fancharts} we note the following:
\begin{itemize}
\item Whilst for models CBD, RH, M7 and PLAT the fans at age 85 are wider than at age 65 in accordance with historical evidence (see \citet[Appendix B]{Cairns2011a}), for models LC and APC the fans at age 85 are narrower than at age 65. This suggest that forecasts from models LC and APC are not plausible for the dataset used in this paper.
\item Forecasts for the PLAT model show an implausible increase of mortality rates. This is because the central trend is linked to the estimated cohort effect $\gamma_{t-x}$ for the PLAT model (see Figure \ref{fig:CohortPLAT}) which shows a steep upward trend between 1935 and 1955.
\item The central trend and uncertainty levels produced by the each of the models have noticeably differences. This highlights the importance of recognising model risk as a significant issue when modelling mortality \citep{Cairns2011a}.  
\end{itemize}

\begin{figure}[!tb]  
    \centering
    \subfloat[LC]{
    <<plotLCfan, fig.width=4, fig.height=5,out.width="4.5cm", echo=FALSE, warning=FALSE, cache=TRUE, dependson="simAll">>=
library(fanplot)
probs = c(2.5, 10, 25, 50, 75, 90, 97.5)
qxt <- Dxt/Ext 
par(mar=c(4.5,4,1,1))
#age 65
plot(LCfit$years, qxt["65", ], xlim = c(1960, 2061), ylim = c(0.0025,0.2),
     xlab = "year", ylab = "mortality rate (log scale)", pch = 20, log = "y")
fan(t(LCsim$rates["65", , ]), start = 2012,
    probs = probs, n.fan = 4,
    fan.col = colorRampPalette(c("black", "white")), ln = NULL)
#age 75
points(LCfit$years, qxt["75", ], pch = 20)
fan(t(LCsim$rates["75", , ]), start = 2012,
    probs = probs, n.fan = 4,
    fan.col = colorRampPalette(c("red", "white")), ln = NULL)
#age 75
points(LCfit$years, qxt["85",], pch = 20)
fan(t(LCsim$rates["85", , ]), start = 2012,
    probs = probs, n.fan = 4,
    fan.col = colorRampPalette(c("blue", "white")), ln = NULL)
#labels
text(1965, qxt[c("65","75", "85"),"1990"], labels=c("x=65", "x=75", "x=85"))
@ 
\label{fig:FanLC}}
\subfloat[CBD]{
<<plotCBDfan, fig.width=4, fig.height=5,out.width="4.5cm", echo=FALSE, cache=TRUE, dependson="plotLCfan">>=
par(mar=c(4.5,4,1,1))
#age 65
plot(CBDfit$years, qxt["65", ], xlim = c(1960, 2061), ylim = c(0.0025,0.2),
     xlab = "year", ylab = "mortality rate (log scale)", pch = 20, log = "y")
fan(t(CBDsim$rates["65", , ]), start = 2012,
    probs = probs, n.fan = 4,
    fan.col = colorRampPalette(c("black", "white")), ln = NULL)
#age 75
points(CBDfit$years, qxt["75", ], pch = 20)
fan(t(CBDsim$rates["75", , ]), start = 2012,
    probs = probs, n.fan = 4,
    fan.col = colorRampPalette(c("red", "white")), ln = NULL)
#age 75
points(CBDfit$years, qxt["85",], pch = 20)
fan(t(CBDsim$rates["85", , ]), start = 2012,
    probs = probs, n.fan = 4,
    fan.col = colorRampPalette(c("blue", "white")), ln = NULL)
#labels
text(1965, qxt[c("65","75", "85"),"1990"], labels=c("x=65", "x=75", "x=85"))
@ 

\label{fig:FanCBD}}
\subfloat[APC]{
    <<plotAPCfan, fig.width=4, fig.height=5,out.width="4.5cm", echo=FALSE, cache=TRUE, dependson="plotLCfan">>=
par(mar=c(4.5,4,1,1))
#age 65
plot(APCfit$years, qxt["65", ], xlim = c(1960, 2061), ylim = c(0.0025,0.2),
     xlab = "year", ylab = "mortality rate (log scale)", pch = 20, log = "y")
fan(t(APCsim$rates["65", , ]), start = 2012,
    probs = probs, n.fan = 4,
    fan.col = colorRampPalette(c("black", "white")), ln = NULL)
#age 75
points(APCfit$years, qxt["75", ], pch = 20)
fan(t(APCsim$rates["75", , ]), start = 2012,
    probs = probs, n.fan = 4,
    fan.col = colorRampPalette(c("red", "white")), ln = NULL)
#age 75
points(APCfit$years, qxt["85",], pch = 20)
fan(t(APCsim$rates["85", , ]), start = 2012,
    probs = probs, n.fan = 4,
    fan.col = colorRampPalette(c("blue", "white")), ln = NULL)
#labels
text(1965, qxt[c("65","75", "85"),"1990"], labels=c("x=65", "x=75", "x=85"))
@ 
\label{fig:FanAPC}}\,
\subfloat[RH]{
    <<plotRHfan, fig.width=4, fig.height=5,out.width="4.5cm", echo=FALSE, cache=TRUE, dependson="plotLCfan">>=
par(mar=c(4.5,4,1,1))
#age 65
plot(RHfit$years, qxt["65", ], xlim = c(1960, 2061), ylim = c(0.0025,0.2),
     xlab = "year", ylab = "mortality rate (log scale)", pch = 20, log = "y")
fan(t(RHsim$rates["65", , ]), start = 2012,
    probs = probs, n.fan = 4,
    fan.col = colorRampPalette(c("black", "white")), ln = NULL)
#age 75
points(RHfit$years, qxt["75", ], pch = 20)
fan(t(RHsim$rates["75", , ]), start = 2012,
    probs = probs, n.fan = 4,
    fan.col = colorRampPalette(c("red", "white")), ln = NULL)
#age 85
points(RHfit$years, qxt["85",], pch = 20)
fan(t(RHsim$rates["85", , ]), start = 2012,
    probs = probs, n.fan = 4,
    fan.col = colorRampPalette(c("blue", "white")), ln = NULL)
#labels
text(1965, qxt[c("65","75", "85"),"1990"], labels=c("x=65", "x=75", "x=85"))
@ 
\label{fig:FanRH}}
\subfloat[M7]{
    <<plotM7fan, fig.width=4, fig.height=5,out.width="4.5cm", echo=FALSE, cache=TRUE, dependson="plotLCfan">>=
par(mar=c(4.5,4,1,1))
#age 65
plot(M7fit$years, qxt["65", ], xlim = c(1960, 2061), ylim = c(0.0025,0.2),
     xlab = "year", ylab = "mortality rate (log scale)", pch = 20, log = "y")
fan(t(M7sim$rates["65", , ]), start = 2012,
    probs = probs, n.fan = 4,
    fan.col = colorRampPalette(c("black", "white")), ln = NULL)
#age 75
points(M7fit$years, qxt["75", ], pch = 20)
fan(t(M7sim$rates["75", , ]), start = 2012,
    probs = probs, n.fan = 4,
    fan.col = colorRampPalette(c("red", "white")), ln = NULL)
#age 85
points(M7fit$years, qxt["85",], pch = 20)
fan(t(M7sim$rates["85", , ]), start = 2012,
    probs = probs, n.fan = 4,
    fan.col = colorRampPalette(c("blue", "white")), ln = NULL)
#labels
text(1965, qxt[c("65","75", "85"),"1990"], labels=c("x=65", "x=75", "x=85"))
@ 
\label{fig:FanM7}}
\subfloat[PLAT]{
    <<plotPLATfan, fig.width=4, fig.height=5,out.width="4.5cm", echo=FALSE, cache=TRUE, dependson="plotLCfan">>=
par(mar=c(4.5,4,1,1))
#age 65
plot(PLATfit$years, qxt["65", ], xlim = c(1960, 2061), ylim = c(0.0025,0.2),
     xlab = "year", ylab = "mortality rate (log scale)", pch = 20, log = "y")
fan(t(PLATsim$rates["65", , ]), start = 2012,
    probs = probs, n.fan = 4,
    fan.col = colorRampPalette(c("black", "white")), ln = NULL)
#age 75
points(PLATfit$years, qxt["75", ], pch = 20)
fan(t(PLATsim$rates["75", , ]), start = 2012,
    probs = probs, n.fan = 4,
    fan.col = colorRampPalette(c("red", "white")), ln = NULL)
#age 85
points(PLATfit$years, qxt["85",], pch = 20)
fan(t(PLATsim$rates["85", , ]), start = 2012,
    probs = probs, n.fan = 4,
    fan.col = colorRampPalette(c("blue", "white")), ln = NULL)
#labels
text(1965, qxt[c("65","75", "85"),"1990"], labels=c("x=65", "x=75", "x=85"))
@ 
\label{fig:FanPLAT}}
		\caption{Fan charts for mortality rates $q_{xt}$ at ages $x=65$ (bottom fan),  $x=75$ (middle fan) and $x=85$ (top fan) from the six models fitted to the England and Wales males population for ages 55-89 and the period 1961-2011. The dots show historical mortality rates for 1961-2011. Shades in the fan represent prediction intervals at the 50\%, 80\% and 95\% level.}\label{fig:Fancharts}		
\end{figure}

\section{Parameter uncertainty and bootstrapping}\label{sec:ParameterUncertainty}

When analysing the uncertainty in mortality projections in an actuarial context it is important to consider all sources of risk. However, the prediction intervals (fan charts) obtained in the previous section only account for the uncertainty arising from the error in the forecast of the period and cohort indexes and ignore the uncertainty arising from the estimation of the parameters of the GAPC model. \\ 

Due to the analytical intractability of many stochastic mortality models, parameter uncertainty is typically accounted for using the bootstrap procedure. This procedure yields $B$ samples  $\alpha_{x}^b, \beta_x^{(1),b},\ldots, \beta_x^{(N),b}, \kappa_t^{(1),b},\ldots,\kappa_t^{(N),b}, \beta_x^{(0),b}, \gamma_{t-x}^b$, $b=1,\ldots,B$, of the parameters of the GAPC model which can then be used to produce confidence and prediction intervals of demographic and actuarial quantities. To sample the parameters of GAPC models we consider the semiparametric bootstrap proposed by \citet{Brouhns2005} and the residual bootstrap first considered in \citet{Koissi2006}:

\begin{itemize}
\item \textbf{Semiparametric bootstrap.} \citet{Brouhns2005} propose a semiparametric bootstrap where first $B$ samples of the number of deaths $d_{xt}^b$, $b=1,\ldots,B$, are generated by sampling from the Poisson Distribution with mean $d_{xt}$. Each bootstrapped sample $d_{xt}^b$, $b=1,\ldots,B$, is then used to re-estimate the model to obtain $B$ bootstrapped parameter estimates $\alpha_{x}^b, \beta_x^{(1),b},\ldots, \beta_x^{(N),b}, \kappa_t^{(1),b},\ldots,\kappa_t^{(N),b}, \beta_x^{(0),b}, \gamma_{t-x}^b$, $b=1,\ldots,B$. \citet{Renshaw2008} use the fitted number of death $\hat{d}_{xt}$ (instead of the observed deaths $d_{xt}$) to perform the sampling from the Poisson distribution. \citet{Wang2005} consider a similar semiparametric approach using a Binomial distribution of deaths.

\item \textbf{Residual bootstrap.} Another possibility is to bootstrap the residuals of the model as suggested by \citet{Koissi2006}. Under this approach the deviance residuals $r_{xt}$ are resampled with replacement to generate $B$ replications $r_{xt}^b$, $b=1,\ldots,B$, which are mapped to the corresponding resampled death counts $d_{xt}^b$, $b=1,\ldots,B$, using the inverse formula.  The model is then re-fitted using these resampled number of deaths to produce $B$ sets of estimated parameters $\alpha_{x}^b, \beta_x^{(1),b},\ldots, \beta_x^{(N),b}, \kappa_t^{(1),b},\ldots,\kappa_t^{(N),b}, \beta_x^{(0),b}, \gamma_{t-x}^b$, $b=1,\ldots,B$. We refer to \citet{Koissi2006} and \citet{Renshaw2008} for details on the inverse formula under a Poisson distribution of deaths and to \citet{Debon2010} for the inverse formula under a Binomial framework. Finally, we note that in implementing the residual bootstrap we follow \citet{Renshaw2008} and apply the inverse formula to produce samples of the observed number of deaths rather than samples of the fitted number of deaths as originally done by \citet{Koissi2006}.  
\end{itemize}

In what follows we illustrate the assessment of \vlad{the} parameter uncertainty using \vlad{the} package \rpkg{StMoMo}. In doing so we deviate from the England and Wales example we have used so far and use instead New Zealand mortality data. This new example follows closely the work of \citet{Li2014} who uses New Zealand mortality data to compare several simulation strategies for assessing the risk in mortality projections with a Poisson Lee-Carter model. The main rationale for the change of dataset is that parameter uncertainty is particularly important when analysing the mortality of smaller populations such as smaller countries or pension plans\footnote{\revI{While the population of England and Wales in 2008 was 54.8 million, the population of New Zealand in 2008 was 4.3 million.}}. This new example also serves as a means for illustrating the use of \vlad{the \rpkg{StMoMo} package} with other datasets. \\

%Reinstate if re-running New Zeland example
% <<HMDpassword, eval=TRUE, echo=FALSE>>=  
% username <- ""
% password <- ""
% @ 

Mortality data for New Zealand can be extracted from the \cite{HumanMortalityDatabase2014} using function \rcode{hmd.mx} of  \vlad{the} \rpkg{demography} package with the code:

%Reinstate if re-running New Zeland example
% <<loadNZ, eval=TRUE, results='hide', cache=TRUE, warning=FALSE, message=FALSE >>=  
%   library(demography)
%   NZdata <- hmd.mx(country = "NZL_NP", username = username,
%                    password = password)
% @ 

<<loadNZ, eval=FALSE>>=  
  library(demography)
  NZdata <- hmd.mx(country = "NZL_NP", username = username,
                   password = password)
@ 

\noindent
We note that the \rcode{username} and \rcode{password} above are for the Human Mortality Database and should be replaced appropriately. Following \citet{Li2014}, we fit a Poisson Lee-Carter model to New Zealand male data for ages 0 to 89 and for the period 1985 to 2008. This can be carried out with the commands:
%Reinstate if re-running New Zeland example
% <<fitNZ, eval=TRUE, results='hide', cache=TRUE, warning=FALSE, message=FALSE, dependson="loadNZ" >>=  
%   Ext_NZ <- NZdata$pop$male 
%   Dxt_NZ <- NZdata$rate$male * Ext_NZ
%   LCfit_NZ <- fit(lc(), Dxt = Dxt_NZ, Ext = Ext_NZ, ages = NZdata$age, 
%              years = NZdata$year, ages.fit = 0:89, years.fit = 1985:2008)
% @ 

<<fitNZ, eval=FALSE>>=  
  Ext_NZ <- NZdata$pop$male 
  Dxt_NZ <- NZdata$rate$male * Ext_NZ
  LCfit_NZ <- fit(lc(), Dxt = Dxt_NZ, Ext = Ext_NZ, ages = NZdata$age, 
             years = NZdata$year, ages.fit = 0:89, years.fit = 1985:2008)
@

In \vlad{the \rpkg{StMoMo} package} the bootstrap of GAPC stochastic mortality models is implemented with the generic function \rcode{bootstrap}. This functions supports both the semiparametric bootstrap and the residual bootstrap. For instance, 5000 semiparametric bootstrap samples of the Lee-Carter model can be obtained with the code\footnote{We note that the bootstrap is a computer intensive procedure. \revI{In particular, the 5000 semiparametric bootstrap samples of the Lee-Carter model took about two hours to run in a computer with an Intel Core i5-3320m processor
running at 2.60 GHz under Windows 7 Home Premium Edition (64 bits) with 8 GB of RAM.}}:
  <<bootLCNZ, eval=FALSE>>=  
LCboot_NZ <- bootstrap(LCfit_NZ, nBoot = 5000, type = "semiparametric")
@

%Reinstate if re-running New Zeland example
%Do this this way as bootstrapping takes a while
% <<loadLCNZboot, eval=TRUE,echo = FALSE, message=FALSE, warning=FALSE, cache=TRUE>>=
% load("LCbootSemiObsNZ.RData")
% LCboot_NZ <- LCbootSemiObs
% @
\noindent
The output from function \rcode{bootstrap} is an object of \vlad{the} class \rcode{"bootStMoMo"} in which the component \rcode{bootParameters} contains the \rcode{nBoot} replications of the bootstrap parameters. A fan chart depicting the 50\%, 80\% and 95\% confidence intervals of the bootstrapped Lee-Carter model (Figure \ref{fig:LCboot}) can be obtained with the command: 
  <<PlotbootLCNZshow, eval=FALSE>>=  
plot(LCboot_NZ)
@
\noindent 
From Figure \ref{fig:LCboot} we note that whilst the parameter uncertainty in the static age function $\alpha_x$ and the period index $\kappa_t^{(1)}$ is modest, the uncertainty in the age-modulating parameters $\beta_x^{(1)}$ is more significant.\\  

\begin{figure}[!tb]  
\centering
%Reinstate if re-running New Zeland example
% <<plotLCNZboot, fig.width=8, fig.height = 2.8, out.width="16cm", echo=FALSE, cache=TRUE, dependson="loadLCNZboot">>=
% plot(LCboot_NZ, nCol = 3)
% @
\includegraphics[width=16cm]{plotLCNZboot.pdf}  
\caption{Bootstrapped parameters for the Poisson Lee-Carter model fitted to the New Zealand male population for ages 0-89 and the period 1985-2008. Shades in the fan represent confidence intervals at the 50\%, 80\% and 95\% level.}\label{fig:LCboot}  
\end{figure}

 Once a stochastic mortality model has been bootstrapped we can simulate it forward to obtain simulated trajectories which account for both the forecast error in the period and cohort indexes and the error in the model fitting. In \rpkg{StMoMo} we can accomplish this using \vlad{the} function \rcode{simulate} applied to an object of class \rcode{"bootStMoMo"}. Thus, to obtain 5000 simulated trajectories of the Lee-Carter model \revI{for the next 24 years taking} into account parameter uncertainty we use the instruction:
% %Reinstate if re-running New Zeland example
% % <<simPULCNZ, cache=TRUE, dependson=c("fitNZ, loadLCNZboot")>>=
% % LCsimPU_NZ <- simulate(LCboot_NZ,  h = 24)
% % @ 
% 
<<simPULCNZ, eval=FALSE>>=
LCsimPU_NZ <- simulate(LCboot_NZ,  h = 24)
@ 
\noindent 
To highlight the impact of parameter risk on mortality rate projections it is instructive to compare  prediction intervals with and without the allowance of parameter uncertainty. \revI{24-year ahead} central forecast together with 5000 trajectories of the Lee-Carter model allowing only for forecast error in the random walk with drift and ignoring the model fitting error can be obtained with the code:
% %Reinstate if re-running New Zeland example
% % <<simLCNZ, cache=TRUE, dependson=c("fitNZ, loadLCNZboot")>>=
% % LCfor_NZ <- forecast(LCfit_NZ, h = 24)
% % LCsim_NZ <- simulate(LCfit_NZ, nsim = 5000, h = 24)
% % @ 
% 
<<simLCNZ, eval=FALSE>>=
LCfor_NZ <- forecast(LCfit_NZ, h = 24)
LCsim_NZ <- simulate(LCfit_NZ, nsim = 5000, h = 24)
@ 
\noindent

Figure \ref{fig:LCNZmxtPred} depicts 95\% prediction intervals for mortality rates at age 40, 60 and 80 with and without allowance for parameter uncertainty. This graph was produced with the code:
<<plotLCPredshow, eval=FALSE>>=  
#Observed, fitted and central forecasts
mxt <- LCfit_NZ$Dxt / LCfit_NZ$Ext
mxtHat <- fitted(LCfit_NZ, type = "rates")           
mxtCentral <- LCfor_NZ$rates
#95% Prediction intervals without parameter uncertainty
mxtPred2.5 <- apply(LCsim_NZ$rates, c(1, 2), quantile, probs = 0.025)
mxtPred97.5 <- apply(LCsim_NZ$rates, c(1, 2), quantile, probs = 0.975)
#95% intervals with parameter uncertainty (in sample, and predictions)
mxtHatPU2.5 <- apply(LCsimPU_NZ$fitted, c(1, 2), quantile, probs = 0.025)  
mxtHatPU97.5 <- apply(LCsimPU_NZ$fitted, c(1, 2), quantile, probs = 0.975)  
mxtPredPU2.5 <- apply(LCsimPU_NZ$rates, c(1, 2), quantile, probs = 0.025)
mxtPredPU97.5 <- apply(LCsimPU_NZ$rates, c(1, 2), quantile, probs = 0.975)
#Plot
x <- c("40", "60", "80")
matplot(LCfit_NZ$years, t(mxt[x, ]), 
        xlim = range(LCfit_NZ$years, LCfor_NZ$years),
        ylim = range(mxtHatPU97.5[x, ], mxtPredPU2.5[x, ], mxt[x, ]),
        type = "p", xlab = "years", ylab = "mortality rates (log scale)",
        log = "y", pch = 20, col = "black")
matlines(LCfit_NZ$years, t(mxtHat[x, ]), lty = 1, col = "black")
matlines(LCfit_NZ$years, t(mxtHatPU2.5[x, ]), lty = 5, col = "red")
matlines(LCfit_NZ$years, t(mxtHatPU97.5[x, ]), lty = 5, col = "red")
matlines(LCfor_NZ$years, t(mxtCentral[x, ]), lty = 4, col = "black")
matlines(LCsim_NZ$years, t(mxtPred2.5[x, ]), lty = 3, col = "black")
matlines(LCsim_NZ$years, t(mxtPred97.5[x, ]), lty = 3, col = "black")
matlines(LCsimPU_NZ$years, t(mxtPredPU2.5[x, ]), lty = 5, col = "red")
matlines(LCsimPU_NZ$years, t(mxtPredPU97.5[x, ]), lty = 5, col = "red")
text(1986, mxtHatPU2.5[x, "1995"], labels = c("x=40", "x=60", "x=80"))
@

\begin{figure}[!tb]  
\centering
%Reinstate if re-running New Zeland example
% <<plotLCPred, eval=TRUE,  fig.width=5, fig.height=5,out.width="7.5cm", results='hide', warning=FALSE, message=FALSE, cache=TRUE, echo=FALSE, dependson=c("simLCNZ", "simPULCNZ", "fitNZ")>>=
% par(mar=c(4.5,4,1,1))
% mxt <- LCfit_NZ$Dxt / LCfit_NZ$Ext
% mxtHat <- fitted(LCfit_NZ, type = "rates")           
% mxtCentral <- LCfor_NZ$rates
% # 95% Prediction intervals without PU
% mxtPred2.5 <- apply(LCsim_NZ$rates, c(1, 2), quantile, probs = 0.025)
% mxtPred97.5 <- apply(LCsim_NZ$rates, c(1, 2), quantile, probs = 0.975)
% # 95% intervals with PU (in sample, and predictions)
% mxtHatPU2.5 <- apply(LCsimPU_NZ$fitted, c(1, 2), quantile, probs = 0.025)  
% mxtHatPU97.5 <- apply(LCsimPU_NZ$fitted, c(1, 2), quantile, probs = 0.975)  
% mxtPredPU2.5 <- apply(LCsimPU_NZ$rates, c(1, 2), quantile, probs = 0.025)
% mxtPredPU97.5 <- apply(LCsimPU_NZ$rates, c(1, 2), quantile, probs = 0.975)
% 
% #Plot
% x <- c("40", "60", "80")
% matplot(LCfit_NZ$years, t(mxt[x, ]), xlim = range(LCfit_NZ$years, LCfor_NZ$years),
%      ylim = range(mxtHatPU97.5[x, ], mxtPredPU2.5[x, ], mxt[x, ]),
%      type = "p", xlab = "years", ylab = "mortality rates (log scale)", 
%      log = "y", pch = 20, col = "black", cex = 0.75)
% matlines(LCfit_NZ$years, t(mxtHat[x, ]), lty = 1, col = "black")
% matlines(LCfit_NZ$years, t(mxtHatPU2.5[x, ]), lty = 5, col ="red")
% matlines(LCfit_NZ$years, t(mxtHatPU97.5[x, ]), lty = 5, col ="red")
% matlines(LCfor_NZ$years, t(mxtCentral[x, ]), lty = 4, col = "black")
% matlines(LCsim_NZ$years, t(mxtPred2.5[x, ]), lty = 3, col = "black", lwd = 2)
% matlines(LCsim_NZ$years, t(mxtPred97.5[x, ]), lty = 3, col = "black", lwd = 2)
% matlines(LCsimPU_NZ$years, t(mxtPredPU2.5[x, ]), lty = 5, col = "red")
% matlines(LCsimPU_NZ$years, t(mxtPredPU97.5[x, ]), lty = 5, col = "red")
% text(1986, mxtHatPU2.5[x, "1995"],
%      labels = c("x = 40",  "x = 60", "x = 80"))
% @ 
\includegraphics[width=7.5cm]{plotLCPred.pdf}    
\caption{95\% Prediction intervals for mortality rates $q_{xt}$ at ages $x=40$ (bottom lines), $x=60$ (middle lines) and $x=80$ (top lines) for the Poisson Lee-Carter model fitted to the New Zealand male population for ages 0-89 and the period 1985-2008. Dots show historical mortality rates for 1985-2008 and solid black lines show the corresponding fitted rates.  Dashed lines represent central forecast, black dotted lines represent 95\% prediction intervals excluding parameter uncertainty and dot-dashed red lines depict 95\% confidence and prediction intervals including parameter uncertainty.}\label{fig:LCNZmxtPred}  
\end{figure}

\noindent
In Figure \ref{fig:LCNZmxtPred} we can clearly see that parameter uncertainty has an important impact on the prediction intervals. This is particularly notable at age 40 where in year 2030, for instance, the width of the prediction interval with parameter uncertainty is around 3 times bigger than that without parameter uncertainty. These results are in line with the results obtained by \citet{Li2014} using the same dataset.\footnote{We note that our prediction intervals \pietro{without and with} parameter uncertainty correspond to  methods (I) and (IVa) in \citet{Li2014}, respectively.}



\section{Conclusion}\label{sec:Conclusions}
In this paper we have introduced the family of Generalised Age-Period-Cohort stochastic mortality models by paralleling the standard framework of generalised linear models.  In addition, we have presented the \R package \rpkg{StMoMo}  which takes advantage of the unifying framework of the GAPC family to provide tools for fitting a diverse number of stochastic mortality models, assessing their goodness of fit and \vlad{also} performing mortality projections. A key feature of the GAPC family and of \rpkg{StMoMo} is that they not only encompass models from the Lee-Carter and CBD families, but can also accommodate possible new models. Furthermore, model risk is a prevalent issue when forecasting mortality and we \vlad{therefore} believe that the possibility of easily implementing and comparing a wide range of models makes \rpkg{StMoMo} a valuable addition to the toolkit for measuring and managing longevity risk.\\

%Furthermore, model risk is a prevalent issue when forcasting mortality and we thus believe that \rpkg{StMoMo}, as a unifying framework  that allows the easy impementation and comparison of a wide range of models, is a valuable addition to the toolkit for measuring and managing longevity risk.\\

Our package can be expanded in several directions. The current version of \rpkg{StMoMo} only allows the use of a log link with a Poisson distribution of deaths or a logit link with a Binomial distribution of deaths. However, we plan to expand  the possible combinations of error distribution and link function to include, for instance, Binomial errors with a complementary log-log link as suggested by \cite{Currie2014}. \\

In agreement with the standard in the literature, \vlad{the \rpkg{StMoMo} package} assumes a multivariate random walk with drift for the dynamics of the period index $\boldsymbol{\kappa}_t$. Nevertheless, the package could be expanded to accommodate more general time series models which might be more appropriate for some datasets or some models. For instance, \citet{Renshaw2006} use an ARIMA$(2,1,0)$ for $\kappa_t^{(1)}$ when forecasting mortality for England and Wales males with the Lee-Carter model with cohorts, while \citet{Plat2009a} uses general univariate ARIMA$(p,d,q)$ for each $\kappa_t^{(i)}$, $i = 1,\ldots, N$, fitted using seemingly unrelated regressions to allow for correlations between the different period indexes. \\

In the GAPC family it is assumed that the age-modulating terms $\beta^{(i)}_x$, $i=0,\ldots, N$, are either non-parametric functions of age which need to be estimated or parametric functions of age with a pre-specified functional form $f^{(i)}(x)$. This latter case could be extended to include the more general case of a pre-specified functional form with a set of parameters that need to be estimated, that is, $\beta^{(i)}_x = f^{(i)}(x;\theta_i)$ with $\theta_i$ being some model parameters. This generalisation would allow the implementation of the family of models considered in \citet{Hunt2014}.\\

Finally, the increasing attention that multiple population mortality models are receiving in the mortality forecasting literature provides an avenue for a wealth of extensions of our package. In particular, we are \revI{considering the development of} a two-population version of \rpkg{StMoMo} that implements the class of relative two-population models considered in \citet{Villegas2015}.


\section{Acknowledgements}\label{sec:Acknowledgements}
This work stems from research funded by the Institute \& Faculty of Actuaries and the Life \& Longevity Markets Association. We thank Sveinn Gunnlaugsson and Mario Sison for their assistance in the development and testing of the preliminary code underlying \rpkg{StMoMo}. We are indebted to Steven Haberman for his encouragement in developing this project and for commenting on a previous version of this paper.  We also thank Anastasios Bardoutsos for his suggestions on the structure of the package and Andrew Hunt for many useful discussions which have helped improved this work. Andr\'es Villegas thanks Ana Deb\'on for first showing him the great things that can be done with the \rpkg{R} package \rpkg{gnm} and selflessly sharing her code and knowledge.


 \bibliographystyle{elsarticle-harv} % outcomment this 
% \bibliography{library} % if more than one, comma separated
\begin{thebibliography}{47}
\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi
\expandafter\ifx\csname url\endcsname\relax
  \def\url#1{\texttt{#1}}\fi
\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL }\fi

\bibitem[{Alai and Sherris(2014)}]{Alai2014b}
Alai, D.~H., Sherris, M., 2014. {Rethinking age-period-cohort mortality trend
  models}. Scandinavian Actuarial Journal~(3), 208--227.

\bibitem[{Aro and Pennanen(2011)}]{Aro2011}
Aro, H., Pennanen, T., 2011. {A user-friendly approach to stochastic mortality
  modelling}. European Actuarial Journal 1, 151--167.

\bibitem[{Booth et~al.(2014)Booth, Hyndman, and Tickle}]{Booth2014}
Booth, H., Hyndman, R.~J., Tickle, L., 2014. {Prospective life tables}. In:
  Charpentier, A. (Ed.), Computational Actuarial Science with R. CRC Press,
  Boca Raton, FL, pp. 319--344.

\bibitem[{Booth et~al.(2002)Booth, Maindonald, and Smith}]{Booth2002}
Booth, H., Maindonald, J., Smith, L., 2002. {Applying Lee-Carter under
  conditions of variable mortality decline.} Population Studies 56~(3),
  325--36.

\bibitem[{B\"{o}rger et~al.(2013)B\"{o}rger, Fleischer, and
  Kuksin}]{Borger2013}
B\"{o}rger, M., Fleischer, D., Kuksin, N., 2013. {Modeling the mortality trend
  under modern solvency regimes}. ASTIN Bulletin 44~(1), 1--38.

\bibitem[{Brouhns et~al.(2005)Brouhns, Denuit, and {Van
  Keilegom}}]{Brouhns2005}
Brouhns, N., Denuit, M., {Van Keilegom}, I., 2005. {Bootstrapping the Poisson
  log-bilinear model for mortality forecasting}. Scandinavian Actuarial
  Journal~(3), 212--224.

\bibitem[{Brouhns et~al.(2002)Brouhns, Denuit, and Vermunt}]{Brouhns2002}
Brouhns, N., Denuit, M., Vermunt, J., 2002. {A Poisson log-bilinear regression
  approach to the construction of projected lifetables}. Insurance: Mathematics
  and Economics 31~(3), 373--393.

\bibitem[{Butt et~al.(2014)Butt, Haberman, and Shang}]{Butt2014}
Butt, Z., Haberman, S., Shang, H.~L., 2014. {ilc: Lee-Carter Mortality Models
  using Iterative Fitting Algorithms}.
\newline\urlprefix\url{http://cran.r-project.org/package=ilc}

\bibitem[{Cairns et~al.(2006)Cairns, Blake, and Dowd}]{Cairns2006a}
Cairns, A.~J.~G., Blake, D., Dowd, K., 2006. {A two-factor model for stochastic
  mortality with parameter uncertainty: theory and calibration}. Journal of
  Risk and Insurance 73~(4), 687--718.

\bibitem[{Cairns et~al.(2011)Cairns, Blake, Dowd, Coughlan, Epstein, and
  Khalaf-Allah}]{Cairns2011a}
Cairns, A.~J.~G., Blake, D., Dowd, K., Coughlan, G.~D., Epstein, D., Khalaf-Allah,
  M., 2011. {Mortality density forecasts: An analysis of six stochastic
  mortality models}. Insurance: Mathematics and Economics 48~(3), 355--367.

\bibitem[{Cairns et~al.(2009)Cairns, Blake, Dowd, Coughlan, Epstein, Ong, and
  Balevich}]{Cairns2009}
Cairns, A.~J.~G., Blake, D., Dowd, K., Coughlan, G.~D., Epstein, D., Ong, A.,
  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.

\bibitem[{Camarda(2012)}]{Camarda2012}
Camarda, C., 2012. {MortalitySmooth: An R package for smoothing Poisson counts
  with P-splines}. Journal of Statistical Software 50~(1), 1--24.

\bibitem[{Clayton and Schifflers(1987)}]{Clayton1987}
Clayton, D., Schifflers, E., 1987. {Models for temporal variation in cancer
  rates. II: Age-period-cohort models.} Statistics in Medicine 6~(4), 469--81.

\bibitem[{Currie(2006)}]{Currie2006}
Currie, I.~D., 2006. {Smoothing and forecasting mortality rates with
  P-splines}.
\newline\urlprefix\url{http://www.macs.hw.ac.uk/~iain/research/talks/Mortality.pdf}

\bibitem[{Currie(2014)}]{Currie2014}
Currie, I.~D., 2014. {On fitting generalized linear and non-linear models of
  mortality}. Scandinavian Actuarial Journal.

\bibitem[{Deb\'{o}n et~al.(2010)Deb\'{o}n, Mart\'{\i}nez-Ruiz, and
  Montes}]{Debon2010}
Deb\'{o}n, A., Mart\'{\i}nez-Ruiz, F., Montes, F., 2010. {A geostatistical
  approach for dynamic life tables: The effect of mortality on remaining
  lifetime and annuities}. Insurance: Mathematics and Economics 47~(3),
  327--336.

\bibitem[{Haberman et~al.(2014)Haberman, Kaishev, Millosaovich, Villegas,
  Baxter, Gaches, Gunnlaugsson, and Sison}]{Haberman2014}
Haberman, S., Kaishev, V.~K., Millosaovich, P., Villegas, A.~M., Baxter, S.,
  Gaches, A., Gunnlaugsson, S., Sison, M., 2014. {Longevity Basis Risk: A
  methodology for assessing basis risk}. Institute and Faculty of Actuaries
  Sessional Research Paper.
\newline\urlprefix\url{http://www.actuaries.org.uk/sites/all/files/documents/pdf/ifoa-llma-longevity-basis-risk-report.pdf}

\bibitem[{Haberman and Renshaw(2009)}]{Haberman2009}
Haberman, S., Renshaw, A., 2009. {On age-period-cohort parametric mortality
  rate projections}. Insurance: Mathematics and Economics 45~(2), 255--270.

\bibitem[{Haberman and Renshaw(2011)}]{Haberman2011a}
Haberman, S., Renshaw, A., 2011. {A comparative study of parametric mortality
  projection models}. Insurance: Mathematics and Economics 48~(1), 35--55.

\bibitem[{Hobcraft et~al.(1982)Hobcraft, Menken, and Preston}]{Hobcraft1982}
Hobcraft, J., Menken, J., Preston, S., 1982. {Age, period, and cohort effects
  in demography: a review.} Population Index 48~(1), 4--43.

\bibitem[{{Human Mortality Database}(2014)}]{HumanMortalityDatabase2014}
{Human Mortality Database}, 2014. {University of California, Berkeley (USA),
  and Max Planck Institute for Demographic Research (Germany)}.
\newline\urlprefix\url{www.mortality.org}

\bibitem[{Hunt and Blake(2014)}]{Hunt2014}
Hunt, A., Blake, D., 2014. {A general procedure for constructing mortality
  models}. North American Actuarial Journal 18~(1), 116--138.

\bibitem[{Hunt and Blake(2015)}]{Hunt2015b}
Hunt, A., Blake, D., 2015. {On the Structure and Classification of Mortality
  Models}. Pension Institute Working Paper.
\newline\urlprefix\url{http://www.pensions-institute.org/workingpapers/wp1506.pdf}

\bibitem[{Hunt and Villegas(2015)}]{Hunt2014a}
Hunt, A., Villegas, A.~M., 2015. {Robustness and convergence in the Lee-Carter
  model with cohorts}. Insurance: Mathematics and Economics 64, 186--202.

\bibitem[{Hyndman(2014)}]{Hyndman2014}
Hyndman, R.~J., 2014. {demography: Forecasting mortality, fertility, migration
  and population data}.
\newline\urlprefix\url{http://cran.r-project.org/package=demography}

\bibitem[{Hyndman and Khandakar(2008)}]{Hyndman2008}
Hyndman, R.~J., Khandakar, Y., 2008. {Automatic time series forecasting: the
  forecast package for R}. Journal of Statistical Software 27~(3), 1--22.

\bibitem[{Hyndman and Ullah(2007)}]{Hyndman2007}
Hyndman, R.~J., Ullah, S., 2007. {Robust forecasting of mortality and fertility
  rates: A functional data approach}. Computational Statistics \& Data Analysis
  51~(10), 4942--4956.

\bibitem[{Koissi et~al.(2006)Koissi, Shapiro, and Hognas}]{Koissi2006}
Koissi, M., Shapiro, A., Hognas, G., 2006. {Evaluating and extending the
  Lee-Carter model for mortality forecasting: Bootstrap confidence interval}.
  Insurance: Mathematics and Economics 38~(1), 1--20.

\bibitem[{Lee and Miller(2001)}]{Lee2001}
Lee, R., Miller, T., 2001. {Evaluating the performance of the Lee-Carter method
  for forecasting mortality.} Demography 38~(4), 537--49.

\bibitem[{Lee and Carter(1992)}]{Lee1992}
Lee, R.~D., Carter, L.~R., 1992. {Modeling and forecasting U.S. mortality}.
  Journal of the American Statistical Association 87~(419), 659--671.

\bibitem[{Li(2014)}]{Li2014}
Li, J., 2014. {A quantitative comparison of simulation strategies for mortality
  projection}. Annals of Actuarial Science 8~(2), 281--297.

\bibitem[{Lov\'{a}sz(2011)}]{Lovasz2011}
Lov\'{a}sz, E., 2011. {Analysis of Finnish and Swedish mortality data with
  stochastic mortality models}. European Actuarial Journal 1~(2), 259--289.

\bibitem[{McCullagh and Nelder(1989)}]{McCullagh1989}
McCullagh, P., Nelder, J., 1989. {Generalized Linear Models}, 2nd Edition.
  Chapman \& Hal, London.

\bibitem[{Oeppen and Vaupel(2002)}]{Oeppen2002}
Oeppen, J., Vaupel, J.~W., 2002. {Broken limits to life expectancy}. Science
  296~(5570), 1029--1031.

\bibitem[{O'Hare and Li(2012)}]{OHare2012}
O'Hare, C., Li, Y., 2012. {Explaining young mortality}. Insurance: Mathematics
  and Economics 50~(1), 12--25.

\bibitem[{Plat(2009)}]{Plat2009a}
Plat, R., 2009. {On stochastic mortality modeling}. Insurance: Mathematics and
  Economics 45~(3), 393--404.

\bibitem[{{R Core Team}(2014)}]{RCoreTeam2014}
{R Core Team}, 2014. {R: A Language and Environment for Statistical Computing}.
  R Foundation for Statistical Computing, Vienna, Austria.
\newline\urlprefix\url{http://www.r-project.org/}

\bibitem[{Renshaw and Haberman(2003)}]{Renshaw2003a}
Renshaw, A., Haberman, S., 2003. {Lee-Carter mortality forecasting with
  age-specific enhancement}. Insurance: Mathematics and Economics 33~(2),
  255--272.

\bibitem[{Renshaw and Haberman(2006)}]{Renshaw2006}
Renshaw, A., Haberman, S., 2006. {A cohort-based extension to the Lee-Carter
  model for mortality reduction factors}. Insurance: Mathematics and Economics
  38~(3), 556--570.

\bibitem[{Renshaw and Haberman(2008)}]{Renshaw2008}
Renshaw, A., Haberman, S., 2008. {On simulation-based approaches to risk
  measurement in mortality with specific reference to Poisson Lee-Carter
  modelling}. Insurance: Mathematics and Economics 42~(2), 797--816.

\bibitem[{Renshaw et~al.(1996)Renshaw, Haberman, and Hatzopoulos}]{Renshaw1996}
Renshaw, A., Haberman, S., Hatzopoulos, P., 1996. {The modelling of recent
  mortality trends in United Kingdom male assured lives}. British Actuarial
  Journal 2~(2), 449--477.

\bibitem[{Turner and Firth(2012)}]{Turner2012}
Turner, H., Firth, D., 2012. {Generalized nonlinear models in R: An overview of
  the gnm package}.
\newline\urlprefix\url{http://cran.r-project.org/package=gnm}

\bibitem[{van Berkum et~al.(2014)van Berkum, Antonio, and
  Vellekoop}]{VanBerkum2014}
van Berkum, F., Antonio, K., Vellekoop, M., 2014. {The impact of multiple
  structural changes on mortality predictions}. Scandinavian Actuarial Journal.

\bibitem[{Villegas et~al.(2015)Villegas, Haberman, Kaishev, and
  Millossovich}]{Villegas2015}
Villegas, A.~M., Haberman, S., Kaishev, V., Millossovich, P., 2015. {A
  comparative study of two-population models for the assessment of basis risk
  in longevity hedges}. Working paper.

\bibitem[{Wang and Lu(2005)}]{Wang2005}
Wang, D., Lu, P., 2005. {Modelling and forecasting mortality distributions in
  England and Wales using the Lee-Carter model}. Journal of Applied
  Statistics 32~(9), 873--885.

\bibitem[{Wikipedia(2014)}]{Wikipedia2015}
Wikipedia, 2014. {King Momo, in Wikipedia, The Free Encyclopedia}.
\newline\urlprefix\url{http://en.wikipedia.org/w/index.php?title=King\_Momo\&oldid=635124783}

\bibitem[{Willets(2004)}]{Willets2004}
Willets, R., 2004. {The cohort effect: Insights and explanations}. British
  Actuarial Journal 10~(4), 833--877.

\end{thebibliography}



\end{document}
