| Title: | Stochastic Mortality Modelling |
|---|---|
| Description: | Implementation of the family of generalised age-period-cohort stochastic mortality models. This family of models encompasses many models proposed in the actuarial and demographic literature including the Lee-Carter (1992) <doi:10.2307/2290201> and the Cairns-Blake-Dowd (2006) <doi:10.1111/j.1539-6975.2006.00195.x> models. It includes functions for fitting mortality models, analysing their goodness-of-fit and performing mortality projections and simulations. |
| Authors: | Andres Villegas <[email protected]>, Pietro Millossovich <[email protected]>, Vladimir Kaishev <[email protected]> |
| Maintainer: | Andres Villegas <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 0.4.1.9000 |
| Built: | 2024-11-22 03:54:20 UTC |
| Source: | https://github.com/amvillegas/stmomo |
Utility function to initialise a StMoMo object representing an
Age-Period-Cohort mortality model.
apc(link = c("log", "logit"))apc(link = c("log", "logit"))
link |
defines the link function and random component associated with
the mortality model. |
The created model is either a log-Poisson or a logit-Binomial version of the classical age-period-cohort mortality model which has predictor structure
To ensure identifiability we follow Cairns et al. (2009) and impose constraints
and
An object of class "StMoMo".
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.
APC <- apc() wxt <- genWeightMat(EWMaleData$ages, EWMaleData$years, clip = 3) APCfit <- fit(APC, data = EWMaleData, wxt = wxt) plot(APCfit, parametricbx = FALSE, nCol = 3)APC <- apc() wxt <- genWeightMat(EWMaleData$ages, EWMaleData$years, clip = 3) APCfit <- fit(APC, data = EWMaleData, wxt = wxt) plot(APCfit, parametricbx = FALSE, nCol = 3)
bootstrap is a generic function for bootstrapping Stochastic
Mortality Models. The function invokes particular methods which depend on
the class of the first argument.
bootstrap(object, nBoot, ...)bootstrap(object, nBoot, ...)
object |
an object used to select a method. Typically of class
|
nBoot |
number of bootstrap samples to produce. |
... |
arguments to be passed to or from other methods. |
bootstrap is a generic function which means that new fitting
strategies can be added for particular stochastic mortality models. See
for instance bootstrap.fitStMoMo.
Produce bootstrap parameters of a Stochastic Mortality Model to account for parameter uncertainty.
## S3 method for class 'fitStMoMo' bootstrap(object, nBoot = 1, type = c("semiparametric", "residual"), deathType = c("observed", "fitted"), ...)## S3 method for class 'fitStMoMo' bootstrap(object, nBoot = 1, type = c("semiparametric", "residual"), deathType = c("observed", "fitted"), ...)
object |
an object of class |
nBoot |
number of bootstrap samples to produce. |
type |
type of bootstrapping approach to be applied.
|
deathType |
type of deaths to sample in the semiparametric bootstrap.
|
... |
arguments to be passed to or from other methods. |
When type is "residual" the residual bootstrapping approach
described in Renshaw and Haberman (2008) is applied, which is an
adaptation of the approach of Koissi et al (2006). In the case of a
"logit" link with Binomial responses the adaptation described in
Debon et al, (2010, section 3) is used.
When type is "semiparametric" the semiparametric approach
described in Brouhns et al.(2005) is used. In the case of a "logit"
link with Binomial responses a suitable adaptation is applied. If
deathType is "observed" then the observed deaths are used in
the sampling as in Brouhns et al. (2005) while if deathType is
"fitted" the fitted deaths are used in the sampling as in
Renshaw and Haberman (2008).
A list with class "bootStMoMo" with components:
bootParameters |
a list of of length |
model |
the model fit that has been bootstrapped. |
type |
type of bootstrapping approach applied. |
deathType |
type of deaths sampled in case of semiparametric bootstrap. |
Brouhns, N., Denuit M., & Van Keilegom, I. (2005). Bootstrapping the Poisson log-bilinear model for mortality forecasting. Scandinavian Actuarial Journal, 2005(3), 212-224.
Debon, A., Martinez-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.
Renshaw, A. E., & 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.
simulate.bootStMoMo, plot.bootStMoMo
#Long computing times ## Not run: LCfit <- fit(lc(), data = EWMaleData) LCResBoot <- bootstrap(LCfit, nBoot = 500, type = "residual") plot(LCResBoot) LCSemiObsBoot <- bootstrap(LCfit, nBoot = 500, type = "semiparametric") plot(LCSemiObsBoot) LCSemiFitBoot <- bootstrap(LCfit, nBoot = 500, type = "semiparametric", deathType = "fitted") plot(LCSemiFitBoot) ## End(Not run)#Long computing times ## Not run: LCfit <- fit(lc(), data = EWMaleData) LCResBoot <- bootstrap(LCfit, nBoot = 500, type = "residual") plot(LCResBoot) LCSemiObsBoot <- bootstrap(LCfit, nBoot = 500, type = "semiparametric") plot(LCSemiObsBoot) LCSemiFitBoot <- bootstrap(LCfit, nBoot = 500, type = "semiparametric", deathType = "fitted") plot(LCSemiFitBoot) ## End(Not run)
Utility function to initialise a StMoMo object representing a
Cairns-Blake-Dowd mortality model.
cbd(link = c("logit", "log"))cbd(link = c("logit", "log"))
link |
defines the link function and random component associated with
the mortality model. |
The created model is either a logit-Binomial or a log-Poisson version of the Cairns-Blake-Dowd mortality model which has predictor structure
where is the average age in the data.
An object of class "StMoMo".
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.
StMoMo, central2initial,
m6, m7, m8
CBD <- cbd() CBDfit <- fit(CBD, data = central2initial(EWMaleData), ages.fit = 55:89) plot(CBDfit, parametricbx = FALSE)CBD <- cbd() CBDfit <- fit(CBD, data = central2initial(EWMaleData), ages.fit = 55:89) plot(CBDfit, parametricbx = FALSE)
Transform StMoMoData from central to initial exposures. Initial exposures are computed by adding one half of the deaths to the central exposures.
central2initial(data)central2initial(data)
data |
StMoMoData object of type "central" created with function
|
A StMoMoData object of type "initial".
CBD <- cbd() CBDfit <- fit(CBD, data = central2initial(EWMaleData), ages.fit = 55:89) plot(CBDfit, parametricbx = FALSE)CBD <- cbd() CBDfit <- fit(CBD, data = central2initial(EWMaleData), ages.fit = 55:89) plot(CBDfit, parametricbx = FALSE)
Extract coefficients from a fitted Stochastic Mortality Model
## S3 method for class 'fitStMoMo' coef(object, ...)## S3 method for class 'fitStMoMo' coef(object, ...)
object |
an object of class |
... |
other arguments. |
A list of model parameters with components:
ax |
Vector with the fitted values of the static age function
|
bx |
Matrix with the values of the period age-modulating functions
|
kt |
Matrix with the values of the fitted period indexes
|
b0x |
Vector with the values of the cohort age-modulating function
|
gc |
Vector with the fitted cohort index |
APCfit <- fit(apc(), data = EWMaleData) coef(APCfit)APCfit <- fit(apc(), data = EWMaleData) coef(APCfit)
Age-specific deaths and exposures for England and Wales from the Human Mortality Database. This is an object of class StMoMoData.
EWMaleDataEWMaleData
A list with the following components:
matrix of deaths data.
matrix of exposures data (mid year population estimates).
vector of ages.
vector of years.
the type of exposure in the data (central).
name of the extracted seriesin this case males.
label of the data.
EWMaleData contains deaths and exposures for England and
Wales males for the period 1961-2011 and for ages 0-100.
Data taken from the Human Mortality Database on 5 November 2014.
Human Mortality Database http://www.mortality.org/.
Human Mortality Database (2014). University of California, Berkeley (USA), and Max Planck Institute for Demographic Research (Germany). Available at www.mortality.org.
Extract cohorts from an age-period array. This is useful to
construct a life table or to perform actuarial/demographic
calculations on a cohort basis using the output of several functions
in StMoMo.
extractCohort(A, age = as.numeric(dimnames(A)[[1]][1]), period = as.numeric(dimnames(A)[[2]][1]), cohort = period - age)extractCohort(A, age = as.numeric(dimnames(A)[[1]][1]), period = as.numeric(dimnames(A)[[2]][1]), cohort = period - age)
A |
an age-period array with a demographic quantity. This array can have two or more dimensions, with the first dimension being the age and the second dimension being the period (calendar year). Note that the names of these two dimension are taken to represent the possible ages and periods in the array. |
age |
optional age for defining the cohort to be extracted. If
argument |
period |
optional period (calendar year) for defining the cohort to be
extracted. If argument |
cohort |
optional cohort to be extracted. If this argument is provided
then arguments |
If the the input array is two dimensional the the output is a
a vector with the quantity along the cohort. Otherwise if A is an
N-dimensional array the output is an (N-1)-dimensional array with the first
dimension representing the cohort.
fitted.fitStMoMo, forecast.fitStMoMo,
simulate.fitStMoMo, simulate.bootStMoMo
LCfit <- fit(lc(), data = EWMaleData, ages.fit = 55:89) #Plot forecast mortality rates for the 1950 cohort LCfor <- forecast(LCfit) plot(55:61, extractCohort(fitted(LCfit, type = "rates"), cohort = 1950), type = "l", log = "y", xlab = "age", ylab = "Mortality rate", main = "Mortality rates for the 1950 cohort", xlim = c(55,89), ylim = c(0.005, 0.12)) lines(62:89, extractCohort(LCfor$rates, cohort = 1950), lty = 2, col = "blue") #Plot 10 simulated sets of mortality rates for the cohort # aged 60 in year 2010 (i.e., the 1950 cohort) LCsim <- simulate(LCfit, nsim = 10) mSim <- extractCohort(LCsim$rates, age = 60, period = 2010) plot(55:61, extractCohort(fitted(LCfit, type = "rates"), cohort = 1950), type = "l", log = "y", xlab = "age", ylab = "Mortality rate", main = "Mortality rates for the 1950 cohort", xlim = c(55,89), ylim = c(0.005, 0.12)) matlines(62:89, mSim, lty = 2)LCfit <- fit(lc(), data = EWMaleData, ages.fit = 55:89) #Plot forecast mortality rates for the 1950 cohort LCfor <- forecast(LCfit) plot(55:61, extractCohort(fitted(LCfit, type = "rates"), cohort = 1950), type = "l", log = "y", xlab = "age", ylab = "Mortality rate", main = "Mortality rates for the 1950 cohort", xlim = c(55,89), ylim = c(0.005, 0.12)) lines(62:89, extractCohort(LCfor$rates, cohort = 1950), lty = 2, col = "blue") #Plot 10 simulated sets of mortality rates for the cohort # aged 60 in year 2010 (i.e., the 1950 cohort) LCsim <- simulate(LCfit, nsim = 10) mSim <- extractCohort(LCsim$rates, age = 60, period = 2010) plot(55:61, extractCohort(fitted(LCfit, type = "rates"), cohort = 1950), type = "l", log = "y", xlab = "age", ylab = "Mortality rate", main = "Mortality rates for the 1950 cohort", xlim = c(55,89), ylim = c(0.005, 0.12)) matlines(62:89, mSim, lty = 2)
fit is a generic function for fitting Stochastic Mortality Models.
The function invokes particular methods which depend on the class of the
first argument.
fit(object, ...)fit(object, ...)
object |
an object used to select a method. Typically of class
|
... |
arguments to be passed to or from other methods. |
fit is a generic function which means that new fitting strategies
can be added for particular stochastic mortality models. See for instance
fit.StMoMo.
Fit a Renshaw and Haberman (Lee-Carter with cohorts) mortality model using the iterative Newton-Raphson procedure presented in Algorithm 1 of Hunt and Villegas (2015). This approach helps solve the well-known robustness and converges issues of the Lee-Carter model with cohort-effects.
## S3 method for class 'rh' fit(object, data = NULL, Dxt = NULL, Ext = NULL, ages = NULL, years = NULL, ages.fit = NULL, years.fit = NULL, oxt = NULL, wxt = NULL, start.ax = NULL, start.bx = NULL, start.kt = NULL, start.b0x = NULL, start.gc = NULL, verbose = TRUE, tolerance = 1e-04, iterMax = 10000, ...)## S3 method for class 'rh' fit(object, data = NULL, Dxt = NULL, Ext = NULL, ages = NULL, years = NULL, ages.fit = NULL, years.fit = NULL, oxt = NULL, wxt = NULL, start.ax = NULL, start.bx = NULL, start.kt = NULL, start.b0x = NULL, start.gc = NULL, verbose = TRUE, tolerance = 1e-04, iterMax = 10000, ...)
object |
an object of class |
data |
an optional object of type StMoMoData containing information on
deaths and exposures to be used for fitting the model. This is typically created
with function |
Dxt |
optional matrix of deaths data. |
Ext |
optional matrix of observed exposures of the same dimension of
|
ages |
optional vector of ages corresponding to rows of |
years |
optional vector of years corresponding to rows of |
ages.fit |
optional vector of ages to include in the fit. Must be a
subset of |
years.fit |
optional vector of years to include in the fit. Must be a
subset of |
oxt |
optional matrix/vector or scalar of known offset to be used in fitting the model. This can be used to specify any a priori known component to be added to the predictor during fitting. |
wxt |
optional matrix of 0-1 weights to be used in the fitting process.
This can be used, for instance, to zero weight some cohorts in the data.
See |
start.ax |
optional vector with starting values for |
start.bx |
optional matrix with starting values for |
start.kt |
optional matrix with starting values for |
start.b0x |
optional vector with starting values for |
start.gc |
optional vector with starting values for |
verbose |
a logical value. If |
tolerance |
a positive numeric value specifying the tolerance level for convergence. |
iterMax |
a positive integer specifying the maximum number of iterations to perform. |
... |
arguments to be passed to or from other methods. |
model |
the object of class |
ax |
vector with the fitted values of the static age function
|
bx |
matrix with the values of the period age-modulating functions
|
kt |
matrix with the values of the fitted period indexes
|
b0x |
vector with the values of the cohort age-modulating function
|
gc |
vector with the fitted cohort index |
data |
StMoMoData object provided for fitting the model. |
Dxt |
matrix of deaths used in the fitting. |
Ext |
matrix of exposures used in the fitting. |
oxt |
matrix of known offset values used in the fitting. |
wxt |
matrix of 0-1 weights used in the fitting. |
ages |
vector of ages used in the fitting. |
years |
vector of years used in the fitting. |
cohorts |
vector of cohorts used in the fitting. |
fittingModel |
output from the iterative fitting algorithm. |
loglik |
log-likelihood of the model. If the fitting failed to
converge this is set to |
deviance |
deviance of the model. If the fitting failed to
converge this is set to |
npar |
effective number of parameters in the model. If the fitting
failed to converge this is set to |
nobs |
number of observations in the model fit. If the fitting
failed to converge this is set to |
fail |
|
conv |
|
Hunt, A., & Villegas, A. M. (2015). Robustness and convergence in the Lee-Carter model with cohorts. Insurance: Mathematics and Economics, 64, 186-202.
LCfit <- fit(lc(), data = EWMaleData, ages.fit = 55:89) wxt <- genWeightMat(55:89, EWMaleData$years, clip = 3) RHfit <- fit(rh(), data = EWMaleData, ages.fit = 55:89, wxt = wxt, start.ax = LCfit$ax, start.bx = LCfit$bx, start.kt = LCfit$kt) plot(RHfit) #Impose approximate constraint as in Hunt and Villegas (2015) ## Not run: RHapprox <- rh(approxConst = TRUE) RHapproxfit <- fit(RHapprox, data = EWMaleData, ages.fit = 55:89, wxt = wxt) plot(RHapproxfit) ## End(Not run)LCfit <- fit(lc(), data = EWMaleData, ages.fit = 55:89) wxt <- genWeightMat(55:89, EWMaleData$years, clip = 3) RHfit <- fit(rh(), data = EWMaleData, ages.fit = 55:89, wxt = wxt, start.ax = LCfit$ax, start.bx = LCfit$bx, start.kt = LCfit$kt) plot(RHfit) #Impose approximate constraint as in Hunt and Villegas (2015) ## Not run: RHapprox <- rh(approxConst = TRUE) RHapproxfit <- fit(RHapprox, data = EWMaleData, ages.fit = 55:89, wxt = wxt) plot(RHapproxfit) ## End(Not run)
Fit a Stochastic Mortality Model to a given data set. The fitting is done
using package gnm.
## S3 method for class 'StMoMo' fit(object, data = NULL, Dxt = NULL, Ext = NULL, ages = NULL, years = NULL, ages.fit = NULL, years.fit = NULL, oxt = NULL, wxt = NULL, start.ax = NULL, start.bx = NULL, start.kt = NULL, start.b0x = NULL, start.gc = NULL, verbose = TRUE, ...)## S3 method for class 'StMoMo' fit(object, data = NULL, Dxt = NULL, Ext = NULL, ages = NULL, years = NULL, ages.fit = NULL, years.fit = NULL, oxt = NULL, wxt = NULL, start.ax = NULL, start.bx = NULL, start.kt = NULL, start.b0x = NULL, start.gc = NULL, verbose = TRUE, ...)
object |
an object of class |
data |
an optional object of type StMoMoData containing information on
deaths and exposures to be used for fitting the model. This is typically created
with function |
Dxt |
optional matrix of deaths data. |
Ext |
optional matrix of observed exposures of the same dimension of
|
ages |
optional vector of ages corresponding to rows of |
years |
optional vector of years corresponding to rows of |
ages.fit |
optional vector of ages to include in the fit. Must be a
subset of |
years.fit |
optional vector of years to include in the fit. Must be a
subset of |
oxt |
optional matrix/vector or scalar of known offset to be used in fitting the model. This can be used to specify any a priori known component to be added to the predictor during fitting. |
wxt |
optional matrix of 0-1 weights to be used in the fitting process.
This can be used, for instance, to zero weight some cohorts in the data.
See |
start.ax |
optional vector with starting values for |
start.bx |
optional matrix with starting values for |
start.kt |
optional matrix with starting values for |
start.b0x |
optional vector with starting values for |
start.gc |
optional vector with starting values for |
verbose |
a logical value. If |
... |
arguments to be passed to or from other methods. This can be
used to control the fitting parameters of |
Fitting is done using function gnm within package
gnm. This is equivalent to minimising (maximising) the deviance
(log-likelihood) of the model. Ages and years in the data should be of
type numeric. Data points with zero exposure are assigned a zero weight
and are ignored in the fitting process. Similarly, NA are assigned a
zero weight and ignored in the fitting process. Parameter estimates can be
plotted using function plot.fitStMoMo.
A list with class "fitStMoMo" with components:
model |
the object of class |
ax |
vector with the fitted values of the static age function
|
bx |
matrix with the values of the period age-modulating functions
|
kt |
matrix with the values of the fitted period indexes
|
b0x |
vector with the values of the cohort age-modulating function
|
gc |
vector with the fitted cohort index |
data |
StMoMoData object provided for fitting the model. |
Dxt |
matrix of deaths used in the fitting. |
Ext |
matrix of exposures used in the fitting. |
oxt |
matrix of known offset values used in the fitting. |
wxt |
matrix of 0-1 weights used in the fitting. |
ages |
vector of ages used in the fitting. |
years |
vector of years used in the fitting. |
cohorts |
vector of cohorts used in the fitting. |
fittingModel |
output from the call to |
loglik |
log-likelihood of the model. If the fitting failed to
converge this is set to |
deviance |
deviance of the model. If the fitting failed to
converge this is set to |
npar |
effective number of parameters in the model. If the fitting
failed to converge this is set to |
nobs |
number of observations in the model fit. If the fitting
failed to converge this is set to |
fail |
|
conv |
|
@seealso StMoMoData, genWeightMat,
plot.fitStMoMo, EWMaleData
# Lee-Carter model only for older ages LCfit <- fit(lc(), data = EWMaleData, ages.fit = 55:89) plot(LCfit) # Use arguments Dxt, Ext, ages, years to pass fitting data LCfit <- fit(lc(), Dxt = EWMaleData$Dxt, Ext = EWMaleData$Ext, ages = EWMaleData$ages, years = EWMaleData$years, ages.fit = 55:89) plot(LCfit) # APC model weigthing out the 3 first and last cohorts wxt <- genWeightMat(EWMaleData$ages, EWMaleData$years, clip = 3) APCfit <- fit(apc(), data = EWMaleData, wxt = wxt) plot(APCfit, parametricbx = FALSE, nCol = 3) # Set verbose = FALSE for silent fitting APCfit <- fit(apc(), data = EWMaleData, wxt = wxt, verbose = FALSE) ## Not run: # Poisson Lee-Carter model with the static age function set to # the mean over time of the log-death rates constLCfix_ax <- function(ax, bx, kt, b0x, gc, wxt, ages){ c1 <- sum(bx, na.rm = TRUE) bx <- bx / c1 kt <- kt * c1 list(ax = ax, bx = bx, kt = kt, b0x = b0x, gc = gc) } LCfix_ax <- StMoMo(link = "log", staticAgeFun = FALSE, periodAgeFun = "NP", constFun = constLCfix_ax) LCfix_axfit <- fit(LCfix_ax, data= EWMaleData, oxt = rowMeans(log(EWMaleData$Dxt / EWMaleData$Ext))) plot(LCfix_axfit) ## End(Not run)# Lee-Carter model only for older ages LCfit <- fit(lc(), data = EWMaleData, ages.fit = 55:89) plot(LCfit) # Use arguments Dxt, Ext, ages, years to pass fitting data LCfit <- fit(lc(), Dxt = EWMaleData$Dxt, Ext = EWMaleData$Ext, ages = EWMaleData$ages, years = EWMaleData$years, ages.fit = 55:89) plot(LCfit) # APC model weigthing out the 3 first and last cohorts wxt <- genWeightMat(EWMaleData$ages, EWMaleData$years, clip = 3) APCfit <- fit(apc(), data = EWMaleData, wxt = wxt) plot(APCfit, parametricbx = FALSE, nCol = 3) # Set verbose = FALSE for silent fitting APCfit <- fit(apc(), data = EWMaleData, wxt = wxt, verbose = FALSE) ## Not run: # Poisson Lee-Carter model with the static age function set to # the mean over time of the log-death rates constLCfix_ax <- function(ax, bx, kt, b0x, gc, wxt, ages){ c1 <- sum(bx, na.rm = TRUE) bx <- bx / c1 kt <- kt * c1 list(ax = ax, bx = bx, kt = kt, b0x = b0x, gc = gc) } LCfix_ax <- StMoMo(link = "log", staticAgeFun = FALSE, periodAgeFun = "NP", constFun = constLCfix_ax) LCfix_axfit <- fit(LCfix_ax, data= EWMaleData, oxt = rowMeans(log(EWMaleData$Dxt / EWMaleData$Ext))) plot(LCfix_axfit) ## End(Not run)
Returns fitted values for the data used in fitting a Stochastic Mortality Model.
## S3 method for class 'fitStMoMo' fitted(object, type = c("link", "rates", "deaths"), ...)## S3 method for class 'fitStMoMo' fitted(object, type = c("link", "rates", "deaths"), ...)
object |
an object of class |
type |
the type of the fitted values that should be returned. The
alternatives are |
... |
other arguments. |
A matrix with the fitted values.
LCfit <- fit(lc(), data = EWMaleData, ages.fit = 55:89) matplot(LCfit$ages, fitted(LCfit), type = "l", lty = 1, col = rainbow(length(LCfit$years)), xlab = "year", ylab = "log death rate", main = "Fitted rates") uxthat <- fitted(LCfit, type = "rates") uxt <- LCfit$Dxt / LCfit$Ext plot(LCfit$years, uxt["65", ], xlab = "year", ylab = "death rate", main = "fitted vs. observed rates at age 65") lines(LCfit$years, uxthat["65", ])LCfit <- fit(lc(), data = EWMaleData, ages.fit = 55:89) matplot(LCfit$ages, fitted(LCfit), type = "l", lty = 1, col = rainbow(length(LCfit$years)), xlab = "year", ylab = "log death rate", main = "Fitted rates") uxthat <- fitted(LCfit, type = "rates") uxt <- LCfit$Dxt / LCfit$Ext plot(LCfit$years, uxt["65", ], xlab = "year", ylab = "death rate", main = "fitted vs. observed rates at age 65") lines(LCfit$years, uxthat["65", ])
Forecast mortality rates using a Stochastic Mortality Model fit.
The period indexes are forecasted
using ether a Multivariate Random Walk with Drift (MRWD) or
independent ARIMA models. The cohort index
is forecasted using an ARIMA.
By default an ARIMA with a constant is used.
## S3 method for class 'fitStMoMo' forecast(object, h = 50, level = c(80, 95), oxt = NULL, gc.order = c(1, 1, 0), gc.include.constant = TRUE, jumpchoice = c("fit", "actual"), kt.method = c("mrwd", "iarima"), kt.order = NULL, kt.include.constant = TRUE, kt.lookback = NULL, gc.lookback = NULL, ...)## S3 method for class 'fitStMoMo' forecast(object, h = 50, level = c(80, 95), oxt = NULL, gc.order = c(1, 1, 0), gc.include.constant = TRUE, jumpchoice = c("fit", "actual"), kt.method = c("mrwd", "iarima"), kt.order = NULL, kt.include.constant = TRUE, kt.lookback = NULL, gc.lookback = NULL, ...)
object |
an object of class |
h |
number of years ahead to forecast. |
level |
confidence level for prediction intervals of the period and cohort indices. |
oxt |
optional matrix/vector or scalar of known offset to be added in the forecasting. This can be used to specify any a priori known component to be added to the forecasted predictor. |
gc.order |
a specification of the ARIMA model for the cohort effect:
the three components |
gc.include.constant |
a logical value indicating if the ARIMA model
should include a constant value. The default is |
jumpchoice |
option to select the jump-off rates, i.e. the rates
from the final year of observation, to use in projections of mortality
rates. |
kt.method |
optional forecasting method for the period index.
The alternatives are |
kt.order |
an optional matrix with one row per period index
specifying the ARIMA models: for the ith row (ith period index) the three
components |
kt.include.constant |
an optional vector of logical values
indicating if the ARIMA model for the ith period index should include a
constant value. The default is |
kt.lookback |
optional argument to specify the look-back window to use
in the estimation of the time series model for the period indexes. By
default all the estimated values are used. If
|
gc.lookback |
optional argument to specify the look-back window to use
in the estimation of the ARIMA model for the cohort effect. By
default all the estimated values are used in estimating the ARIMA
model. If |
... |
other arguments for |
If kt.method is "mrwd", fitting and forecasting of
the time series model for the period indexes is done with a
Multivariate Random Walk with Drift using the function
mrwd.
If kt.method is "iarima", fitting and forecasting of
the time series model for the period indexes is done with
independent arima models using the function iarima.
See this latter function for details on input arguments
kt.order and kt.include.constant.
Fitting and forecasting of the ARIMA model for the cohort index
is done with function Arima from package
forecast. See the latter function for further details on
input arguments gc.order and gc.include.constant.
Note that in some cases forecast of the cohort effects may be
needed for a horizon longer than h. This is the case when
in the fitted model the most recent cohorts have been zero weighted.
The forecasted cohorts can be seen in gc.f$cohorts.
A list of class "forStMoMo" with components:
rates |
a matrix with the point forecast of the rates. |
ages |
vector of ages corresponding to the rows of |
years |
vector of years for which a forecast has been produced. This
corresponds to the columns of |
kt.f |
forecasts of period indexes of the model. This is a list with
the |
gc.f |
forecasts of cohort index of the model. This is a list with
the |
oxt.f |
the offset used in the forecast. |
fitted |
a matrix with the fitted in-sample rates of the model for the years for which the mortality model was fitted. |
model |
the model fit from which the forecast was produced. |
jumpchoice |
Jump-off method used in the forecast. |
kt.method |
method used in the forecast of the period index. |
#Lee-Carter (random walk with drift) LCfit <- fit(lc(), data = EWMaleData, ages.fit = 55:89) LCfor <- forecast(LCfit) plot(LCfor) #Lee-Carter (forecast with ARIMA(1,1,2) with drift ) LCfor.iarima1 <- forecast(LCfit, kt.method = "iarima", kt.order = c(1, 1 , 2)) plot(LCfor.iarima1) #Lee-Carter (forecast with auto.arima) LCfor.iarima2 <- forecast(LCfit, kt.method = "iarima") plot(LCfor.iarima2) #CBD (Multivariate random walk with drift) CBDfit <- fit(cbd(), data = central2initial(EWMaleData), ages.fit = 55:89) CBDfor <- forecast(CBDfit) plot(CBDfor, parametricbx = FALSE) #CBD (Independent Arima models) kt.order <- matrix(c(1, 1, 2, #ARIMA(1, 1, 2) for k[1] 0, 1, 1), #ARIMA(0, 1, 1) for k[2] nrow = 2, ncol = 3, byrow = TRUE) CBDfor.iarima <- forecast(CBDfit, kt.method = "iarima", kt.order = kt.order) plot(CBDfor.iarima, parametricbx = FALSE) #APC: Compare forecast with different models for the cohort index wxt <- genWeightMat(55:89, EWMaleData$years, clip = 3) APCfit <- fit(apc(), data = EWMaleData, ages.fit = 55:89, wxt = wxt) APCfor1 <- forecast(APCfit) plot(APCfor1, parametricbx = FALSE, nCol = 3) APCfor2 <- forecast(APCfit, gc.order = c(0, 2, 2)) plot(APCfor2, only.gc = TRUE) plot(c(APCfit$years, APCfor1$years), cbind(APCfor1$fitted, APCfor1$rates)["65", ], type = "l", xlab = "year", ylab = "Mortality rate at age 65", main = "Forecasts with different models for gc") lines(APCfor2$years, APCfor2$rates["65", ], col = "blue") points(APCfit$years, (APCfit$Dxt / APCfit$Ext)["65", ], pch = 19) #Compare Lee-Carter forecast using: # 1. Fitted jump-off rates and all history for kt # 2. Actual jump-off rates and all history for kt # 3. Fitted jump-off rates and only history for # the past 30 years of kt (i.e 1982-2011) LCfor1 <- forecast(LCfit) LCfor2 <- forecast(LCfit, jumpchoice = "actual") LCfor3 <- forecast(LCfit, kt.lookback = 30) plot(LCfit$years, (LCfit$Dxt / LCfit$Ext)["60", ], xlim = range(LCfit$years, LCfor1$years), ylim = range((LCfit$Dxt / LCfit$Ext)["60", ], LCfor1$rates["60", ], LCfor2$rates["60", ], LCfor3$rates["60", ]), type = "p", xlab = "year", ylab = "rate", main = "Lee-Carter: Forecast of mortality rates at age 60") lines(LCfit$years, fitted(LCfit, type = "rates")["60", ]) lines(LCfor1$years, LCfor1$rates["60", ], lty = 2) lines(LCfor2$years, LCfor2$rates["60", ], lty = 3, col = "blue") lines(LCfor3$years, LCfor3$rates["60", ], lty = 4, col = "red") legend("topright",legend = c("Fitted jump-off", "Actual jump-off", "Fitted jump-off, 30 year look-back"), lty = 1:3, col = c("black", "blue", "red"))#Lee-Carter (random walk with drift) LCfit <- fit(lc(), data = EWMaleData, ages.fit = 55:89) LCfor <- forecast(LCfit) plot(LCfor) #Lee-Carter (forecast with ARIMA(1,1,2) with drift ) LCfor.iarima1 <- forecast(LCfit, kt.method = "iarima", kt.order = c(1, 1 , 2)) plot(LCfor.iarima1) #Lee-Carter (forecast with auto.arima) LCfor.iarima2 <- forecast(LCfit, kt.method = "iarima") plot(LCfor.iarima2) #CBD (Multivariate random walk with drift) CBDfit <- fit(cbd(), data = central2initial(EWMaleData), ages.fit = 55:89) CBDfor <- forecast(CBDfit) plot(CBDfor, parametricbx = FALSE) #CBD (Independent Arima models) kt.order <- matrix(c(1, 1, 2, #ARIMA(1, 1, 2) for k[1] 0, 1, 1), #ARIMA(0, 1, 1) for k[2] nrow = 2, ncol = 3, byrow = TRUE) CBDfor.iarima <- forecast(CBDfit, kt.method = "iarima", kt.order = kt.order) plot(CBDfor.iarima, parametricbx = FALSE) #APC: Compare forecast with different models for the cohort index wxt <- genWeightMat(55:89, EWMaleData$years, clip = 3) APCfit <- fit(apc(), data = EWMaleData, ages.fit = 55:89, wxt = wxt) APCfor1 <- forecast(APCfit) plot(APCfor1, parametricbx = FALSE, nCol = 3) APCfor2 <- forecast(APCfit, gc.order = c(0, 2, 2)) plot(APCfor2, only.gc = TRUE) plot(c(APCfit$years, APCfor1$years), cbind(APCfor1$fitted, APCfor1$rates)["65", ], type = "l", xlab = "year", ylab = "Mortality rate at age 65", main = "Forecasts with different models for gc") lines(APCfor2$years, APCfor2$rates["65", ], col = "blue") points(APCfit$years, (APCfit$Dxt / APCfit$Ext)["65", ], pch = 19) #Compare Lee-Carter forecast using: # 1. Fitted jump-off rates and all history for kt # 2. Actual jump-off rates and all history for kt # 3. Fitted jump-off rates and only history for # the past 30 years of kt (i.e 1982-2011) LCfor1 <- forecast(LCfit) LCfor2 <- forecast(LCfit, jumpchoice = "actual") LCfor3 <- forecast(LCfit, kt.lookback = 30) plot(LCfit$years, (LCfit$Dxt / LCfit$Ext)["60", ], xlim = range(LCfit$years, LCfor1$years), ylim = range((LCfit$Dxt / LCfit$Ext)["60", ], LCfor1$rates["60", ], LCfor2$rates["60", ], LCfor3$rates["60", ]), type = "p", xlab = "year", ylab = "rate", main = "Lee-Carter: Forecast of mortality rates at age 60") lines(LCfit$years, fitted(LCfit, type = "rates")["60", ]) lines(LCfor1$years, LCfor1$rates["60", ], lty = 2) lines(LCfor2$years, LCfor2$rates["60", ], lty = 3, col = "blue") lines(LCfor3$years, LCfor3$rates["60", ], lty = 4, col = "red") legend("topright",legend = c("Fitted jump-off", "Actual jump-off", "Fitted jump-off, 30 year look-back"), lty = 1:3, col = c("black", "blue", "red"))
Returns forecasts and other information for a group of independent arima series.
## S3 method for class 'iarima' forecast(object, h = 10, level = c(80, 95), fan = FALSE, ...)## S3 method for class 'iarima' forecast(object, h = 10, level = c(80, 95), fan = FALSE, ...)
object |
an object of class |
h |
Number of periods for forecasting. |
level |
confidence level for prediction intervals. |
fan |
if |
... |
other arguments. |
An object of class "iarimaForecast" with components:
model |
a list containing information about the fitted arima models. |
mean |
array with the central forecast. |
lower |
three dimensional array with lower limits for prediction intervals. |
upper |
three dimensional array with upper limits for prediction intervals. |
level |
the confidence values associated with the prediction intervals. |
@export
Returns forecasts and other information for a Multivariate Random Walk with Drift model.
## S3 method for class 'mrwd' forecast(object, h = 10, level = c(80, 95), fan = FALSE, ...)## S3 method for class 'mrwd' forecast(object, h = 10, level = c(80, 95), fan = FALSE, ...)
object |
an object of class |
h |
Number of periods for forecasting. |
level |
confidence level for prediction intervals. |
fan |
if |
... |
other arguments. |
An object of class "mrwdForecast" with components:
model |
a list containing information about the fitted model. |
mean |
array with the central forecast. |
lower |
three dimensional array with lower limits for prediction intervals. |
upper |
three dimensional array with upper limits for prediction intervals. |
level |
the confidence values associated with the prediction intervals. |
@export
Generates a weight matrix given a group of ages and years
and a set of cohorts which are to be given zero weight. This
is useful for excluding some data points when fitting a
Stochastic Mortality Model (see fit.StMoMo).
genWeightMat(ages, years, clip = 0, zeroCohorts = NULL)genWeightMat(ages, years, clip = 0, zeroCohorts = NULL)
ages |
vector of ages. |
years |
vector of years. |
clip |
number of cohorts in the boundary to assign a zero weight. This can be be used to zero weigh some of the first and last cohorts in the data. |
zeroCohorts |
other cohort for which a zero weight is to be assigned. |
A 0-1 matrix with 0 for the zero-weighed cohorts.
#Zero-weight the first three and last three cohorts wxt1 <- genWeightMat(55:89, EWMaleData$years, clip = 3) APCfit1 <- fit(apc(), data = EWMaleData, ages.fit = 55:89, wxt = wxt1) plot(APCfit1, parametricbx = FALSE, nCol = 3) #Also Zero-weight the 1886 cohort wxt2 <- genWeightMat(55:89, EWMaleData$years, clip = 3, zeroCohorts = 1886) APCfit2 <- fit(apc(), data = EWMaleData, ages.fit = 55:89, wxt = wxt2) plot(APCfit2, parametricbx = FALSE, nCol = 3)#Zero-weight the first three and last three cohorts wxt1 <- genWeightMat(55:89, EWMaleData$years, clip = 3) APCfit1 <- fit(apc(), data = EWMaleData, ages.fit = 55:89, wxt = wxt1) plot(APCfit1, parametricbx = FALSE, nCol = 3) #Also Zero-weight the 1886 cohort wxt2 <- genWeightMat(55:89, EWMaleData$years, clip = 3, zeroCohorts = 1886) APCfit2 <- fit(apc(), data = EWMaleData, ages.fit = 55:89, wxt = wxt2) plot(APCfit2, parametricbx = FALSE, nCol = 3)
Fits independent arima series to x, a multivariate
time series.
iarima(x, order = NULL, include.constant = TRUE, ...)iarima(x, order = NULL, include.constant = TRUE, ...)
x |
numeric matrix with a multivariate time series. Series are arranged in rows with columns representing time. |
order |
an optional matrix with one row per time series
specifying the ARIMA models: for the ith row the three
components |
include.constant |
an optional vector of logical values
indicating if the ARIMA model for the ith series should include a
constant value. The default is |
... |
additional parameters for |
The fitting of the ARIMA models for each time series is done with
function Arima from package
forecast. See the latter function for further details on
input arguments kt.order and kt.include.constant.
an object of class "iarima" with components:
models |
a list with the arima models fitted to each time series. |
x |
the original time series. |
Transform StMoMoData from initial to central exposures. Central exposures are computed by substracting one half of the deaths from the initial exposures.
initial2central(data)initial2central(data)
data |
StMoMoData object of type "initial" created with function
|
A StMoMoData object of type "central".
Utility function to initialise a StMoMo object representing a
Lee-Carter model.
lc(link = c("log", "logit"), const = c("sum", "last", "first"))lc(link = c("log", "logit"), const = c("sum", "last", "first"))
link |
defines the link function and random component associated with
the mortality model. |
const |
defines the constraint to impose to the period index of the
model to ensure identifiability. The alternatives are
|
The created model is either a log-Poisson (see Brouhns et al (2002)) or a logit-Binomial version of the Lee-Carter model which has predictor structure
To ensure identifiability one of the following constraints is imposed
depending on the value of const, and
An object of class "StMoMo".
Brouhns, N., Denuit, M., & Vermunt, J. K. (2002). A Poisson log-bilinear regression approach to the construction of projected lifetables. Insurance: Mathematics and Economics, 31(3), 373-393.
Lee, R. D., & Carter, L. R. (1992). Modeling and forecasting U.S. mortality. Journal of the American Statistical Association, 87(419), 659-671.
#sum(kt) = 0 and log link LC1 <- lc() LCfit1<-fit(LC1, data = EWMaleData, ages.fit = 55:89) plot(LCfit1) #kt[1] = 0 and log link LC2 <- lc(const = "first") LCfit2<-fit(LC2, data = EWMaleData, ages.fit = 55:89) plot(LCfit2) #kt[n] = 0 and logit link LC3 <- lc("logit", "last") LCfit3<-fit(LC3, data = EWMaleData, ages.fit = 55:89) plot(LCfit3)#sum(kt) = 0 and log link LC1 <- lc() LCfit1<-fit(LC1, data = EWMaleData, ages.fit = 55:89) plot(LCfit1) #kt[1] = 0 and log link LC2 <- lc(const = "first") LCfit2<-fit(LC2, data = EWMaleData, ages.fit = 55:89) plot(LCfit2) #kt[n] = 0 and logit link LC3 <- lc("logit", "last") LCfit3<-fit(LC3, data = EWMaleData, ages.fit = 55:89) plot(LCfit3)
Returns the log-likelihood of a fitted Stochastic Mortality Model.
## S3 method for class 'fitStMoMo' logLik(object, ...)## S3 method for class 'fitStMoMo' logLik(object, ...)
object |
an object of class |
... |
other arguments. |
The log-likelihood of the fitted model.
Utility function to initialise a StMoMo object representing the
M6 (CBD with cohorts) extension of the Cairns-Blake-Dowd mortality model
introduced in Cairns et al (2009).
m6(link = c("logit", "log"))m6(link = c("logit", "log"))
link |
defines the link function and random component associated with
the mortality model. |
The created model is either a logit-Binomial or a log-Poisson version of the M6 model which has predictor structure
where is the average age in the data.
Identifiability of the model is accomplished by applying parameters constraints
which ensure that the cohort effect fluctuates around zero and has no linear trend. These constraints are applied using the strategy discussed in Appendix A of Haberman and Renshaw (2011).
An object of class "StMoMo".
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.
Haberman, S., & Renshaw, A. (2011). A comparative study of parametric mortality projection models. Insurance: Mathematics and Economics, 48(1), 35-55.
StMoMo, central2initial,
cbd, m7, m8
M6 <- m6() wxt <- genWeightMat(55:89, EWMaleData$years, clip = 3) M6fit <- fit(M6, data = central2initial(EWMaleData), ages.fit = 55:89) plot(M6fit, parametricbx = FALSE)M6 <- m6() wxt <- genWeightMat(55:89, EWMaleData$years, clip = 3) M6fit <- fit(M6, data = central2initial(EWMaleData), ages.fit = 55:89) plot(M6fit, parametricbx = FALSE)
Utility function to initialise a StMoMo object representing the
M7 extension of the Cairns-Blake-Dowd mortality model introduced
in Cairns et al (2009).
m7(link = c("logit", "log"))m7(link = c("logit", "log"))
link |
defines the link function and random component associated with
the mortality model. |
The created model is either a logit-Binomial or a log-Poisson version of the M7 model which has predictor structure
where is the average age in the data and
is the average value of .
Identifiability of the model is accomplished by applying parameters constraints
which ensure that the cohort effect fluctuates around zero and has no linear or quadratic trend. These constraints are applied using the strategy discussed in Appendix A of Haberman and Renshaw (2011).
An object of class "StMoMo".
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.
Haberman, S., & Renshaw, A. (2011). A comparative study of parametric mortality projection models. Insurance: Mathematics and Economics, 48(1), 35-55.
StMoMo, central2initial,
cbd, m6, m8
M7 <- m7() wxt <- genWeightMat(55:89, EWMaleData$years, clip = 3) M7fit <- fit(M7, data = central2initial(EWMaleData), ages.fit = 55:89) plot(M7fit, parametricbx = FALSE)M7 <- m7() wxt <- genWeightMat(55:89, EWMaleData$years, clip = 3) M7fit <- fit(M7, data = central2initial(EWMaleData), ages.fit = 55:89) plot(M7fit, parametricbx = FALSE)
Utility function to initialise a StMoMo object representing the
M8 extension of the Cairns-Blake-Dowd mortality model introduced
in Cairns et al (2009).
m8(link = c("logit", "log"), xc)m8(link = c("logit", "log"), xc)
link |
defines the link function and random component associated with
the mortality model. |
xc |
constant defining the cohort age-modulating parameter. |
The created model is either a logit-Binomial or a log-Poisson version of the M8 model which has predictor structure
where is the average age in the data and is a
predefined constant.
Identifiability of the model is accomplished by applying parameters
constraint
An object of class "StMoMo".
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.
StMoMo, central2initial,
cbd, m6, m7
M8 <- m8(xc = 89) wxt <- genWeightMat(55:89, EWMaleData$years, clip = 3) M8fit <- fit(M8, data = central2initial(EWMaleData), ages.fit = 55:89) plot(M8fit, parametricbx = FALSE)M8 <- m8(xc = 89) wxt <- genWeightMat(55:89, EWMaleData$years, clip = 3) M8fit <- fit(M8, data = central2initial(EWMaleData), ages.fit = 55:89) plot(M8fit, parametricbx = FALSE)
Fits a Multivariate Random Walk with Drift to x, a
multivariate time series.
mrwd(x)mrwd(x)
x |
numeric matrix with a multivariate time series. Series are arranged in rows with columns representing time. |
For further information on the Multivariate Random Walk with drift see Appendix B in Haberman and Renshaw (2011).
an object of class "mrwd" with components:
drift |
a vector with the estimated drift. |
sigma |
a matrix with the estimated variance covariance matrix. |
fitted |
fitted values. |
residuals |
residuals from the fitted model. That is observed minus fitted values. |
x |
the original time series. |
Haberman, S., & Renshaw, A. (2011). A comparative study of parametric mortality projection models. Insurance: Mathematics and Economics, 48(1), 35-55.
Plot fancharts of bootstrapped parameters of a Stochastic Mortality Model
stored in an object of class "bootStMoMo".
## S3 method for class 'bootStMoMo' plot(x, nCol = 2, parametricbx = TRUE, colour = rgb(0, 0, 0), probs = c(2.5, 10, 25, 50, 75, 90, 97.5), ...)## S3 method for class 'bootStMoMo' plot(x, nCol = 2, parametricbx = TRUE, colour = rgb(0, 0, 0), probs = c(2.5, 10, 25, 50, 75, 90, 97.5), ...)
x |
an object of class |
nCol |
number of columns to use in the plot. |
parametricbx |
if |
colour |
colour to use in the fans. |
probs |
probabilities related to percentiles to plot in the fan chart.
The default |
... |
other arguments. |
#Long computing times ## Not run: CBDfit <- fit(cbd(), data = central2initial(EWMaleData), ages.fit = 55:89) CBDResBoot <- bootstrap(CBDfit, nBoot = 500) plot(CBDResBoot) plot(CBDResBoot, parametricbx = FALSE, probs = seq(2.5, 97.5, 2.5)) ## End(Not run)#Long computing times ## Not run: CBDfit <- fit(cbd(), data = central2initial(EWMaleData), ages.fit = 55:89) CBDResBoot <- bootstrap(CBDfit, nBoot = 500) plot(CBDResBoot) plot(CBDResBoot, parametricbx = FALSE, probs = seq(2.5, 97.5, 2.5)) ## End(Not run)
Plot fitted parameters of a stochastic mortality model of class
"fitStMoMo".
## S3 method for class 'fitStMoMo' plot(x, nCol = 2, parametricbx = TRUE, type = "l", ...)## S3 method for class 'fitStMoMo' plot(x, nCol = 2, parametricbx = TRUE, type = "l", ...)
x |
an object of class |
nCol |
number of columns to use in the plot. |
parametricbx |
if |
type |
what type of plot should be drawn. See
|
... |
additional arguments to control graphical appearance.
See |
#Fit and plot a Lee-Carter model LCfit <- fit(lc(), data = EWMaleData, ages.fit = 55:89) plot(LCfit) plot(LCfit, type = "p", pch = 19) #Fit and plot a CBD model CBDfit <- fit(cbd(), data = central2initial(EWMaleData), ages.fit = 55:89) plot(CBDfit) plot(CBDfit, parametricbx = FALSE) plot(CBDfit, nCol = 1, parametricbx = FALSE, lwd = 2)#Fit and plot a Lee-Carter model LCfit <- fit(lc(), data = EWMaleData, ages.fit = 55:89) plot(LCfit) plot(LCfit, type = "p", pch = 19) #Fit and plot a CBD model CBDfit <- fit(cbd(), data = central2initial(EWMaleData), ages.fit = 55:89) plot(CBDfit) plot(CBDfit, parametricbx = FALSE) plot(CBDfit, nCol = 1, parametricbx = FALSE, lwd = 2)
Plot a forecasted Stochastic Mortality Model of class "forStMoMo".
## S3 method for class 'forStMoMo' plot(x, nCol = 2, parametricbx = TRUE, only.kt = FALSE, only.gc = FALSE, colour = "grey60", ...)## S3 method for class 'forStMoMo' plot(x, nCol = 2, parametricbx = TRUE, only.kt = FALSE, only.gc = FALSE, colour = "grey60", ...)
x |
an object of class |
nCol |
number of columns to use in the plot. |
parametricbx |
if |
only.kt |
If |
only.gc |
If |
colour |
colour to use in the prediction intervals. |
... |
additional arguments to control graphical appearance.
See |
wxt <- genWeightMat(55:89, EWMaleData$years, clip = 3) APCfit <- fit(apc(), data = EWMaleData, ages.fit = 55:89, wxt = wxt) APCfor <- forecast(APCfit) plot(APCfor) plot(APCfor, parametricbx = FALSE, nCol = 3) plot(APCfor, only.kt = TRUE) plot(APCfor, only.gc = TRUE, lwd = 2)wxt <- genWeightMat(55:89, EWMaleData$years, clip = 3) APCfit <- fit(apc(), data = EWMaleData, ages.fit = 55:89, wxt = wxt) APCfor <- forecast(APCfit) plot(APCfor) plot(APCfor, parametricbx = FALSE, nCol = 3) plot(APCfor, only.kt = TRUE) plot(APCfor, only.gc = TRUE, lwd = 2)
Plots the deviance residuals of a Stochastic Mortality Model which are
of class "resStMoMo". Three types of plots
are available: scatter plot of residuals by age, period and cohort,
colour map (heatmap) of the residuals, and a black and white signplot
of the residuals.
## S3 method for class 'resStMoMo' plot(x, type = c("scatter", "colourmap", "signplot"), reslim = NULL, plotAge = TRUE, plotYear = TRUE, plotCohort = TRUE, pch = 20, col = NULL, ...)## S3 method for class 'resStMoMo' plot(x, type = c("scatter", "colourmap", "signplot"), reslim = NULL, plotAge = TRUE, plotYear = TRUE, plotCohort = TRUE, pch = 20, col = NULL, ...)
x |
an object of class |
type |
the type of the plot. The alternatives are
|
reslim |
optional numeric vector of length 2, giving the range of the residuals. |
plotAge |
logical value indicating if the age scatter plot should be
produced. This is only used when |
plotYear |
logical value indicating if the calendar year scatter plot
should be produced. This is only used when |
plotCohort |
logical value indicating if the cohort scatter plot
should be produced. This is only used when |
pch |
optional symbol to use for the points in a scatterplot.
This is only used when |
col |
optional colours to use in plotting. If
|
... |
other plotting parameters to be passed to the plotting functions. This can be used to control the appearance of the plots. |
When type = "scatter" scatter plots of the residuals against age,
calendar year and cohort (year of birth) are produced.
When type = "colourmap" a two dimensional colour map of the
residuals is plotted. This is produced using function
image.plot. See image.plot
for further parameters that can be passed to this type of plots.
When type = "signplot" a two dimensional black and white map of the
residuals is plotted with dark grey representing negative residuals and
light grey representing positive residuals. This is produced using
function image.default.
@seealso residuals.fitStMoMo
CBDfit <- fit(cbd(), data = central2initial(EWMaleData), ages.fit = 55:89) CBDres <- residuals(CBDfit) plot(CBDres) plot(CBDres, type = "signplot") plot(CBDres, type = "colourmap")CBDfit <- fit(cbd(), data = central2initial(EWMaleData), ages.fit = 55:89) CBDres <- residuals(CBDfit) plot(CBDres) plot(CBDres, type = "signplot") plot(CBDres, type = "colourmap")
Obtain predictions from a Stochastic Mortality Model fit.
## S3 method for class 'fitStMoMo' predict(object, years, kt = NULL, gc = NULL, oxt = NULL, type = c("link", "rates"), ...)## S3 method for class 'fitStMoMo' predict(object, years, kt = NULL, gc = NULL, oxt = NULL, type = c("link", "rates"), ...)
object |
an object of class |
years |
vector of years for which a prediction is required. |
kt |
matrix of values of the period indexes to use for the prediction.
If the model has any age-period term this argument needs to be provided and
the number of rows in |
gc |
vector of values of the cohort indexes to use for the prediction.
If the model has a cohort effect this argument needs to be provided.
In this case the length of |
oxt |
optional matrix/vector or scalar of known offset to be used in the prediction. |
type |
the type of the predicted values that should be returned. The
alternatives are |
... |
other arguments. |
This function evaluates
for a fitted Stochastic Mortality model.
In producing a prediction the static age function, , and the
age-modulating parameters, , are taken from
the fitted model in object while the period indexes,
, and cohort index, ,
are taken from the function arguments.
This function can be useful, for instance, in producing forecasts of
mortality rates using time series models different to those available
in forecast.fitStMoMo (See examples below).
A matrix with the predicted values.
## Not run: ##M6 Forecast using VARIMA(1,1) model library(MTS) # fit m6 years <- EWMaleData$years ages.fit <- 55:89 M6 <- m6(link = "log") M6fit <- fit(M6, data = EWMaleData, ages.fit = ages.fit) # Forecast kt using VARIMA(1,1) model from MTS h <- 50 kt.M6 <- t(M6fit$kt) kt.M6.diff <- apply(kt.M6, 2, diff) fit.kt.M6.11 <- VARMA(kt.M6.diff, p = 1, q = 1) pred.ktdiff.M6.11 <- VARMApred(fit.kt.M6.11, h = h) pred.kt.M6.11 <- apply(rbind(tail(kt.M6, n = 1), pred.ktdiff.M6.11$pred), 2, cumsum)[-1, ] # set row names years.forecast <- seq(tail(years, 1) + 1, length.out = h) rownames(pred.kt.M6.11) <- years.forecast # plot kt1 plot(x = c(years, years.forecast), y = c(kt.M6[, 1], pred.kt.M6.11[, 1]), col = rep(c("black", "red"), times = c(length(years), h)), xlab = "time", ylab = "k1") plot(x = c(years, years.forecast), y = c(kt.M6[, 2], pred.kt.M6.11[, 2]), col = rep(c("black", "red"), times = c(length(years), h)), xlab = "time", ylab = "k2") # forecast cohort effect # the following cohorts are required: # from 2012 - 89 = 1923 # to 2061 - 55 = 2006 pred.gc.M6 <- forecast(auto.arima(M6fit$gc, max.d = 1), h = h) # use predict to get rates pred.qxt.M6.11 <- predict(object = M6fit, years = years.forecast, kt = t(pred.kt.M6.11), gc = c(tail(M6fit$gc, 34), pred.gc.M6$mean), type = "rates") qxthatM6 <- fitted(M6fit, type = "rates") # plot mortality profile at age 60, 70 and 80 matplot(1961 : 2061, t(cbind(qxthatM6, pred.qxt.M6.11)[c("60", "70", "80"), ]), type = "l", col = "black", xlab = "years", ylab = "rates", lwd = 1.5) ## End(Not run)## Not run: ##M6 Forecast using VARIMA(1,1) model library(MTS) # fit m6 years <- EWMaleData$years ages.fit <- 55:89 M6 <- m6(link = "log") M6fit <- fit(M6, data = EWMaleData, ages.fit = ages.fit) # Forecast kt using VARIMA(1,1) model from MTS h <- 50 kt.M6 <- t(M6fit$kt) kt.M6.diff <- apply(kt.M6, 2, diff) fit.kt.M6.11 <- VARMA(kt.M6.diff, p = 1, q = 1) pred.ktdiff.M6.11 <- VARMApred(fit.kt.M6.11, h = h) pred.kt.M6.11 <- apply(rbind(tail(kt.M6, n = 1), pred.ktdiff.M6.11$pred), 2, cumsum)[-1, ] # set row names years.forecast <- seq(tail(years, 1) + 1, length.out = h) rownames(pred.kt.M6.11) <- years.forecast # plot kt1 plot(x = c(years, years.forecast), y = c(kt.M6[, 1], pred.kt.M6.11[, 1]), col = rep(c("black", "red"), times = c(length(years), h)), xlab = "time", ylab = "k1") plot(x = c(years, years.forecast), y = c(kt.M6[, 2], pred.kt.M6.11[, 2]), col = rep(c("black", "red"), times = c(length(years), h)), xlab = "time", ylab = "k2") # forecast cohort effect # the following cohorts are required: # from 2012 - 89 = 1923 # to 2061 - 55 = 2006 pred.gc.M6 <- forecast(auto.arima(M6fit$gc, max.d = 1), h = h) # use predict to get rates pred.qxt.M6.11 <- predict(object = M6fit, years = years.forecast, kt = t(pred.kt.M6.11), gc = c(tail(M6fit$gc, 34), pred.gc.M6$mean), type = "rates") qxthatM6 <- fitted(M6fit, type = "rates") # plot mortality profile at age 60, 70 and 80 matplot(1961 : 2061, t(cbind(qxthatM6, pred.qxt.M6.11)[c("60", "70", "80"), ]), type = "l", col = "black", xlab = "years", ylab = "rates", lwd = 1.5) ## End(Not run)
Compute deviance residuals of a fitted Stochastic Mortality Model.
These residuals can be plotted using plot.resStMoMo.
## S3 method for class 'fitStMoMo' residuals(object, scale = TRUE, ...)## S3 method for class 'fitStMoMo' residuals(object, scale = TRUE, ...)
object |
an object of class |
scale |
logical indicating whether the residuals should be scaled
or not by dividing the deviance by the overdispersion of the model.
Default is |
... |
other arguments. |
An object of class "resStMoMo" with the residuals. This
object has components:
residuals |
a matrix with the residuals. |
ages |
ages corresponding to the rows in |
years |
years corresponding to the columns in |
CBDfit <- fit(cbd(), data = central2initial(EWMaleData), ages.fit = 55:89) CBDres <- residuals(CBDfit) plot(CBDres)CBDfit <- fit(cbd(), data = central2initial(EWMaleData), ages.fit = 55:89) CBDres <- residuals(CBDfit) plot(CBDres)
Utility function to initialise a StMoMo object representing a
Renshaw and Haberman (Lee-Carter with cohorts) mortality model introduced
in Renshaw and Haberman (2006).
rh(link = c("log", "logit"), cohortAgeFun = c("1", "NP"), approxConst = FALSE)rh(link = c("log", "logit"), cohortAgeFun = c("1", "NP"), approxConst = FALSE)
link |
defines the link function and random component associated with
the mortality model. |
cohortAgeFun |
defines the cohort age modulating parameter
|
approxConst |
defines if the approximate identifiability constraint of
Hunt and Villegas (2015) is applied or not. If |
The created model is either a log-Poisson or a logit-Binomial version of the Renshaw and Haberman model which has predictor structure
or
depending on the value of argument cohortAgeFun.
To ensure identifiability the following constraints are imposed
plus
if cohortAgeFun = "NP"
In addition, if approxConst=TRUE then the approximate
identifiability constraint
is applied to improve the stability and robustness of the model (see Hunt and Villegas (2015)).
By default as this model has shown to be more
stable (see Haberman and Renshaw (2011) and Hunt and Villegas (2015)).
An object of class "StMoMo" or "rh".
Haberman, S., & Renshaw, A. (2011). A comparative study of parametric mortality projection models. Insurance: Mathematics and Economics, 48(1), 35-55.
Hunt, A., & Villegas, A. M. (2015). Robustness and convergence in the Lee-Carter model with cohorts. Insurance: Mathematics and Economics, 64, 186-202.
Renshaw, A. E., & Haberman, S. (2006). A cohort-based extension to the Lee-Carter model for mortality reduction factors. Insurance: Mathematics and Economics, 38(3), 556-570.
LCfit <- fit(lc(), data = EWMaleData, ages.fit = 55:89) wxt <- genWeightMat(55:89, EWMaleData$years, clip = 3) RHfit <- fit(rh(), data = EWMaleData, ages.fit = 55:89, wxt = wxt, start.ax = LCfit$ax, start.bx = LCfit$bx, start.kt = LCfit$kt) plot(RHfit) #Impose approximate constraint as in Hunt and Villegas (2015) ## Not run: RHapprox <- rh(approxConst = TRUE) RHapproxfit <- fit(RHapprox, data = EWMaleData, ages.fit = 55:89, wxt = wxt) plot(RHapproxfit) ## End(Not run)LCfit <- fit(lc(), data = EWMaleData, ages.fit = 55:89) wxt <- genWeightMat(55:89, EWMaleData$years, clip = 3) RHfit <- fit(rh(), data = EWMaleData, ages.fit = 55:89, wxt = wxt, start.ax = LCfit$ax, start.bx = LCfit$bx, start.kt = LCfit$kt) plot(RHfit) #Impose approximate constraint as in Hunt and Villegas (2015) ## Not run: RHapprox <- rh(approxConst = TRUE) RHapproxfit <- fit(RHapprox, data = EWMaleData, ages.fit = 55:89, wxt = wxt) plot(RHapproxfit) ## End(Not run)
Simulate future sample paths from a Bootstrapped Stochastic Mortality Model.
The period indexes are modelled
using ether a Multivariate Random Walk with Drift (MRWD) or
independent ARIMA models. The cohort index
is modelled using an ARIMA.
By default an ARIMA with a constant is used.
## S3 method for class 'bootStMoMo' simulate(object, nsim = 1, seed = NULL, h = 50, oxt = NULL, gc.order = c(1, 1, 0), gc.include.constant = TRUE, jumpchoice = c("fit", "actual"), kt.method = c("mrwd", "iarima"), kt.order = NULL, kt.include.constant = TRUE, kt.lookback = NULL, gc.lookback = NULL, ...)## S3 method for class 'bootStMoMo' simulate(object, nsim = 1, seed = NULL, h = 50, oxt = NULL, gc.order = c(1, 1, 0), gc.include.constant = TRUE, jumpchoice = c("fit", "actual"), kt.method = c("mrwd", "iarima"), kt.order = NULL, kt.include.constant = TRUE, kt.lookback = NULL, gc.lookback = NULL, ...)
object |
an object of class |
nsim |
number of sample paths to simulate from each bootstrapped
sample. Thus if there are |
seed |
either |
h |
number of years ahead to forecast. |
oxt |
optional array/matrix/vector or scalar of known offset to be added in the simulations. This can be used to specify any a priori known component to be added to the simulated predictor. |
gc.order |
a specification of the ARIMA model for the cohort effect:
the three components |
gc.include.constant |
a logical value indicating if the ARIMA model
should include a constant value. The default is |
jumpchoice |
option to select the jump-off rates, i.e. the rates
from the final year of observation, to use in projections of mortality
rates. |
kt.method |
optional forecasting method for the period index.
The alternatives are |
kt.order |
an optional matrix with one row per period index
specifying the ARIMA models: for the ith row (ith period index) the three
components |
kt.include.constant |
an optional vector of logical values
indicating if the ARIMA model for the ith period index should include a
constant value. The default is |
kt.lookback |
optional argument to specify the look-back window to use
in the estimation of the time series model for the period indexes. By
default all the estimated values are used. If
|
gc.lookback |
optional argument to specify the look-back window to use
in the estimation of the ARIMA model for the cohort effect. By
default all the estimated values are used in estimating the ARIMA
model. If |
... |
other arguments. |
For further details see simulate.fitStMoMo.
A list of class "simStMoMo" with components
rates |
a three dimensional array with the future simulated rates. |
ages |
vector of ages corresponding to the first dimension of
|
years |
vector of years for which a simulations has been produced.
This corresponds to the second dimension of |
kt.s |
information on the simulated paths of the period indices of
the model. This is a list with the simulated paths of |
gc.s |
information on the simulated paths of the cohort index of the
model. This is a list with the simulated paths of |
oxt.s |
a three dimensional array with the offset used in the simulations. |
fitted |
a three dimensional array with the in-sample rates of the model for the years for which the mortality model was fitted (and bootstrapped). |
jumpchoice |
Jump-off method used in the simulation. |
kt.method |
method used in the modelling of the period index. |
model |
the bootstrapped model from which the simulations were produced. |
bootstrap.fitStMoMo, simulate.fitStMoMo
#Long computing times ## Not run: #Lee-Carter: Compare projection with and without parameter uncertainty library(fanplot) LCfit <- fit(lc(), data = EWMaleData) LCResBoot <- bootstrap(LCfit, nBoot = 500) LCResBootsim <- simulate(LCResBoot) LCsim <- simulate(LCfit, nsim = 500) plot(LCfit$years, log(LCfit$Dxt / LCfit$Ext)["10", ], xlim = range(LCfit$years, LCsim$years), ylim = range(log(LCfit$Dxt / LCfit$Ext)["10", ], log(LCsim$rates["10", , ])), type = "l", xlab = "year", ylab = "log rate", main = "Mortality rate projection at age 10 with and without parameter uncertainty") fan(t(log(LCResBootsim$rates["10", , ])),start = LCResBootsim$years[1], probs = c(2.5, 10, 25, 50, 75, 90, 97.5), n.fan = 4, fan.col = colorRampPalette(c(rgb(0, 0, 1), rgb(1, 1, 1))), ln = NULL) fan(t(log(LCsim$rates["10", 1:(length(LCsim$years) - 3), ])), start = LCsim$years[1], probs = c(2.5, 10, 25, 50, 75, 90, 97.5), n.fan = 4, fan.col = colorRampPalette(c(rgb(1, 0, 0), rgb(1, 1, 1))), ln = NULL) ## End(Not run)#Long computing times ## Not run: #Lee-Carter: Compare projection with and without parameter uncertainty library(fanplot) LCfit <- fit(lc(), data = EWMaleData) LCResBoot <- bootstrap(LCfit, nBoot = 500) LCResBootsim <- simulate(LCResBoot) LCsim <- simulate(LCfit, nsim = 500) plot(LCfit$years, log(LCfit$Dxt / LCfit$Ext)["10", ], xlim = range(LCfit$years, LCsim$years), ylim = range(log(LCfit$Dxt / LCfit$Ext)["10", ], log(LCsim$rates["10", , ])), type = "l", xlab = "year", ylab = "log rate", main = "Mortality rate projection at age 10 with and without parameter uncertainty") fan(t(log(LCResBootsim$rates["10", , ])),start = LCResBootsim$years[1], probs = c(2.5, 10, 25, 50, 75, 90, 97.5), n.fan = 4, fan.col = colorRampPalette(c(rgb(0, 0, 1), rgb(1, 1, 1))), ln = NULL) fan(t(log(LCsim$rates["10", 1:(length(LCsim$years) - 3), ])), start = LCsim$years[1], probs = c(2.5, 10, 25, 50, 75, 90, 97.5), n.fan = 4, fan.col = colorRampPalette(c(rgb(1, 0, 0), rgb(1, 1, 1))), ln = NULL) ## End(Not run)
Simulate future sample paths from a Stochastic Mortality Model.
The period indexes are modelled
using ether a Multivariate Random Walk with Drift (MRWD) or
independent ARIMA models. The cohort index
is modelled using an ARIMA.
By default an ARIMA with a constant is used.
## S3 method for class 'fitStMoMo' simulate(object, nsim = 1000, seed = NULL, h = 50, oxt = NULL, gc.order = c(1, 1, 0), gc.include.constant = TRUE, jumpchoice = c("fit", "actual"), kt.method = c("mrwd", "iarima"), kt.order = NULL, kt.include.constant = TRUE, kt.lookback = NULL, gc.lookback = NULL, ...)## S3 method for class 'fitStMoMo' simulate(object, nsim = 1000, seed = NULL, h = 50, oxt = NULL, gc.order = c(1, 1, 0), gc.include.constant = TRUE, jumpchoice = c("fit", "actual"), kt.method = c("mrwd", "iarima"), kt.order = NULL, kt.include.constant = TRUE, kt.lookback = NULL, gc.lookback = NULL, ...)
object |
an object of class |
nsim |
number of sample paths to simulate. |
seed |
either |
h |
number of years ahead to forecast. |
oxt |
optional matrix/vector or scalar of known offset to be added in the simulations. This can be used to specify any a priori known component to be added to the simulated predictor. |
gc.order |
a specification of the ARIMA model for the cohort effect:
the three components |
gc.include.constant |
a logical value indicating if the ARIMA model
should include a constant value. The default is |
jumpchoice |
option to select the jump-off rates, i.e. the rates
from the final year of observation, to use in projections of mortality
rates. |
kt.method |
optional forecasting method for the period index.
The alternatives are |
kt.order |
an optional matrix with one row per period index
specifying the ARIMA models: for the ith row (ith period index) the three
components |
kt.include.constant |
an optional vector of logical values
indicating if the ARIMA model for the ith period index should include a
constant value. The default is |
kt.lookback |
optional argument to specify the look-back window to use
in the estimation of the time series model for the period indexes. By
default all the estimated values are used. If
|
gc.lookback |
optional argument to specify the look-back window to use
in the estimation of the ARIMA model for the cohort effect. By
default all the estimated values are used in estimating the ARIMA
model. If |
... |
other arguments. |
If kt.method is "mrwd", fitting and simulation of
the time series model for the period indexes is done with a
Multivariate Random Walk with Drift using the function
mrwd.
If kt.method is "iarima", fitting and simulation of
the time series model for the period indexes is done with
independent arima models using the function iarima.
See this latter function for details on input arguments
kt.order and kt.include.constant.
Fitting and simulation of the ARIMA model for the cohort index
is done with function Arima from package
forecast. See the latter function for further details on
input arguments gc.order and gc.include.constant.
Note that in some cases simulations of the
cohort effects may be needed for a horizon longer than h.
This is the case when in the fitted model the most recent cohorts
have been zero weighted. The simulated cohorts can be seen in
gc.s$cohorts.
A list of class "simStMoMo" with components:
rates |
a three dimensional array with the future simulated rates. |
ages |
vector of ages corresponding to the first dimension of
|
years |
vector of years for which a simulations has been produced.
This corresponds to the second dimension of |
kt.s |
information on the simulated paths of the period indexes
of the model. This is a list with the |
gc.s |
information on the simulated paths of the cohort index of
the model. This is a list with the |
oxt.s |
a three dimensional array with the offset used in the simulations. |
fitted |
a three dimensional array with the in-sample rates of the model for the years for which the mortality model was fitted. |
jumpchoice |
Jump-off method used in the simulation. |
kt.method |
method used in the modelling of the period index. |
model |
the model fit from which the simulations were produced. |
#Lee-Carter LCfit <- fit(lc(), data = EWMaleData, ages.fit = 55:89) LCsim.mrwd <- simulate(LCfit, nsim = 100) LCsim.iarima <- simulate(LCfit, nsim = 100, kt.method = "iarima", kt.order = c(1, 1, 2)) par(mfrow=c(2, 2)) plot(LCfit$years, LCfit$kt[1, ], xlim = range(LCfit$years, LCsim.mrwd$kt.s$years), ylim = range(LCfit$kt, LCsim.mrwd$kt.s$sim), type = "l", xlab = "year", ylab = "kt", main = "Lee-Carter: Simulated paths of the period index kt (mrwd)") matlines(LCsim.mrwd$kt.s$years, LCsim.mrwd$kt.s$sim[1, , ], type = "l", lty = 1) plot(LCfit$years, (LCfit$Dxt / LCfit$Ext)["65", ], xlim = range(LCfit$years, LCsim.mrwd$years), ylim = range((LCfit$Dxt / LCfit$Ext)["65", ], LCsim.mrwd$rates["65", , ]), type = "l", xlab = "year", ylab = "rate", main = "Lee-Carter: Simulated mortality rates at age 65") matlines(LCsim.mrwd$years, LCsim.mrwd$rates["65", , ], type = "l", lty = 1) plot(LCfit$years, LCfit$kt[1, ], xlim = range(LCfit$years, LCsim.iarima$kt.s$years), ylim = range(LCfit$kt, LCsim.iarima$kt.s$sim), type = "l", xlab = "year", ylab = "kt", main = "Lee-Carter: Simulated paths of the period index kt (ARIMA(1, 1, 2))") matlines(LCsim.iarima$kt.s$years, LCsim.iarima$kt.s$sim[1, , ], type = "l", lty = 1) plot(LCfit$years, (LCfit$Dxt / LCfit$Ext)["65", ], xlim = range(LCfit$years, LCsim.iarima$years), ylim = range((LCfit$Dxt / LCfit$Ext)["65", ], LCsim.iarima$rates["65", , ]), type = "l", xlab = "year", ylab = "rate", main = "Lee-Carter: Simulated mortality rates at age 65 (ARIMA(1, 1, 2))") matlines(LCsim.iarima$years, LCsim.iarima$rates["65", , ], type = "l", lty = 1) #APC par(mfrow=c(1, 3)) wxt <- genWeightMat(55:89, EWMaleData$years, clip = 3) APCfit <- fit(apc(), data = EWMaleData, ages.fit = 55:89, wxt = wxt) APCsim <- simulate(APCfit, nsim = 100, gc.order = c(1, 1, 0)) plot(APCfit$years, APCfit$kt[1, ], xlim = range(APCfit$years, APCsim$kt.s$years), ylim = range(APCfit$kt, APCsim$kt.s$sim), type = "l", xlab = "year", ylab = "kt", main = "APC: Simulated paths of the period index kt") matlines(APCsim$kt.s$years, APCsim$kt.s$sim[1, , ], type = "l", lty = 1) plot(APCfit$cohorts, APCfit$gc, xlim = range(APCfit$cohorts, APCsim$gc.s$cohorts), ylim = range(APCfit$gc, APCsim$gc.s$sim, na.rm = TRUE), type = "l", xlab = "year", ylab = "kt", main = "APC: Simulated paths of the cohort index (ARIMA(1,1,0))") matlines(APCsim$gc.s$cohorts, APCsim$gc.s$sim, type = "l", lty = 1) plot(APCfit$years, (APCfit$Dxt / APCfit$Ext)["65", ], xlim = range(APCfit$years, APCsim$years), ylim = range((APCfit$Dxt/APCfit$Ext)["65", ], APCsim$rates["65", , ]), type = "l", xlab = "year", ylab = "rate", main = "APC: Simulated of mortality rates at age 65") matlines(APCsim$years, APCsim$rates["65", , ], type = "l", lty = 1) #Compare LC and APC library(fanplot) par(mfrow=c(1, 1)) plot(LCfit$years, (LCfit$Dxt / LCfit$Ext)["65", ], xlim = range(LCfit$years, LCsim.mrwd$years), ylim = range((LCfit$Dxt / LCfit$Ext)["65", ], LCsim.mrwd$rates["65", , ], APCsim$rates["65", , ]), type = "l", xlab = "year", ylab = "rate", main = "Fan chart of mortality rates at age 65 (LC vs. APC)") fan(t(LCsim.mrwd$rates["65", , ]), start = LCsim.mrwd$years[1], probs = c(2.5, 10, 25, 50, 75, 90, 97.5), n.fan = 4, fan.col = colorRampPalette(c(rgb(1, 0, 0), rgb(1, 1, 1))), ln = NULL) fan(t(APCsim$rates["65", 1:(length(APCsim$years) - 3), ]), start = APCsim$years[1], probs = c(2.5, 10, 25, 50, 75, 90, 97.5), n.fan = 4, fan.col = colorRampPalette(c(rgb(0, 0, 1), rgb(1, 1, 1))), ln = NULL)#Lee-Carter LCfit <- fit(lc(), data = EWMaleData, ages.fit = 55:89) LCsim.mrwd <- simulate(LCfit, nsim = 100) LCsim.iarima <- simulate(LCfit, nsim = 100, kt.method = "iarima", kt.order = c(1, 1, 2)) par(mfrow=c(2, 2)) plot(LCfit$years, LCfit$kt[1, ], xlim = range(LCfit$years, LCsim.mrwd$kt.s$years), ylim = range(LCfit$kt, LCsim.mrwd$kt.s$sim), type = "l", xlab = "year", ylab = "kt", main = "Lee-Carter: Simulated paths of the period index kt (mrwd)") matlines(LCsim.mrwd$kt.s$years, LCsim.mrwd$kt.s$sim[1, , ], type = "l", lty = 1) plot(LCfit$years, (LCfit$Dxt / LCfit$Ext)["65", ], xlim = range(LCfit$years, LCsim.mrwd$years), ylim = range((LCfit$Dxt / LCfit$Ext)["65", ], LCsim.mrwd$rates["65", , ]), type = "l", xlab = "year", ylab = "rate", main = "Lee-Carter: Simulated mortality rates at age 65") matlines(LCsim.mrwd$years, LCsim.mrwd$rates["65", , ], type = "l", lty = 1) plot(LCfit$years, LCfit$kt[1, ], xlim = range(LCfit$years, LCsim.iarima$kt.s$years), ylim = range(LCfit$kt, LCsim.iarima$kt.s$sim), type = "l", xlab = "year", ylab = "kt", main = "Lee-Carter: Simulated paths of the period index kt (ARIMA(1, 1, 2))") matlines(LCsim.iarima$kt.s$years, LCsim.iarima$kt.s$sim[1, , ], type = "l", lty = 1) plot(LCfit$years, (LCfit$Dxt / LCfit$Ext)["65", ], xlim = range(LCfit$years, LCsim.iarima$years), ylim = range((LCfit$Dxt / LCfit$Ext)["65", ], LCsim.iarima$rates["65", , ]), type = "l", xlab = "year", ylab = "rate", main = "Lee-Carter: Simulated mortality rates at age 65 (ARIMA(1, 1, 2))") matlines(LCsim.iarima$years, LCsim.iarima$rates["65", , ], type = "l", lty = 1) #APC par(mfrow=c(1, 3)) wxt <- genWeightMat(55:89, EWMaleData$years, clip = 3) APCfit <- fit(apc(), data = EWMaleData, ages.fit = 55:89, wxt = wxt) APCsim <- simulate(APCfit, nsim = 100, gc.order = c(1, 1, 0)) plot(APCfit$years, APCfit$kt[1, ], xlim = range(APCfit$years, APCsim$kt.s$years), ylim = range(APCfit$kt, APCsim$kt.s$sim), type = "l", xlab = "year", ylab = "kt", main = "APC: Simulated paths of the period index kt") matlines(APCsim$kt.s$years, APCsim$kt.s$sim[1, , ], type = "l", lty = 1) plot(APCfit$cohorts, APCfit$gc, xlim = range(APCfit$cohorts, APCsim$gc.s$cohorts), ylim = range(APCfit$gc, APCsim$gc.s$sim, na.rm = TRUE), type = "l", xlab = "year", ylab = "kt", main = "APC: Simulated paths of the cohort index (ARIMA(1,1,0))") matlines(APCsim$gc.s$cohorts, APCsim$gc.s$sim, type = "l", lty = 1) plot(APCfit$years, (APCfit$Dxt / APCfit$Ext)["65", ], xlim = range(APCfit$years, APCsim$years), ylim = range((APCfit$Dxt/APCfit$Ext)["65", ], APCsim$rates["65", , ]), type = "l", xlab = "year", ylab = "rate", main = "APC: Simulated of mortality rates at age 65") matlines(APCsim$years, APCsim$rates["65", , ], type = "l", lty = 1) #Compare LC and APC library(fanplot) par(mfrow=c(1, 1)) plot(LCfit$years, (LCfit$Dxt / LCfit$Ext)["65", ], xlim = range(LCfit$years, LCsim.mrwd$years), ylim = range((LCfit$Dxt / LCfit$Ext)["65", ], LCsim.mrwd$rates["65", , ], APCsim$rates["65", , ]), type = "l", xlab = "year", ylab = "rate", main = "Fan chart of mortality rates at age 65 (LC vs. APC)") fan(t(LCsim.mrwd$rates["65", , ]), start = LCsim.mrwd$years[1], probs = c(2.5, 10, 25, 50, 75, 90, 97.5), n.fan = 4, fan.col = colorRampPalette(c(rgb(1, 0, 0), rgb(1, 1, 1))), ln = NULL) fan(t(APCsim$rates["65", 1:(length(APCsim$years) - 3), ]), start = APCsim$years[1], probs = c(2.5, 10, 25, 50, 75, 90, 97.5), n.fan = 4, fan.col = colorRampPalette(c(rgb(0, 0, 1), rgb(1, 1, 1))), ln = NULL)
Returns one simulated path of the group of independent
arima series in object.
## S3 method for class 'iarima' simulate(object, nsim = 10, seed = NULL, ...)## S3 method for class 'iarima' simulate(object, nsim = 10, seed = NULL, ...)
object |
An object of class |
nsim |
number of periods for the simulated series. |
seed |
either |
... |
other arguments. |
Returns one simulated path of the Multivariate
Random Walk with Drift model in object.
## S3 method for class 'mrwd' simulate(object, nsim = 10, seed = NULL, ...)## S3 method for class 'mrwd' simulate(object, nsim = 10, seed = NULL, ...)
object |
An object of class |
nsim |
number of periods for the simulated series. |
seed |
either |
... |
other arguments. |
Initialises a StMoMo object which represents a Generalised Age-Period-Cohort Stochastic Mortality Model.
StMoMo.
StMoMo(link = c("log", "logit"), staticAgeFun = TRUE, periodAgeFun = "NP", cohortAgeFun = NULL, constFun = function(ax, bx, kt, b0x, gc, wxt, ages) list(ax = ax, bx = bx, kt = kt, b0x = b0x, gc = gc))StMoMo(link = c("log", "logit"), staticAgeFun = TRUE, periodAgeFun = "NP", cohortAgeFun = NULL, constFun = function(ax, bx, kt, b0x, gc, wxt, ages) list(ax = ax, bx = bx, kt = kt, b0x = b0x, gc = gc))
link |
defines the link function and random component associated with
the mortality model. |
staticAgeFun |
logical value indicating if a static age function
|
periodAgeFun |
a list of length |
cohortAgeFun |
defines the cohort age modulating parameter
|
constFun |
function defining the identifiability constraints of the
model. It must be a function of the form
|
R implementation of the family of Generalised Age-Period-Cohort stochastic mortality models. This family of models encompasses many models proposed in the literature including the well-known Lee-Carter model, CBD model and APC model.
StMoMo defines an abstract representation of a Generalised
Age-Period-Cohort (GAPC) Stochastic model that fits within the
general class of generalised non-linear models defined as follows
where
is a static age function;
, are age/period terms;
is the age/cohort term; and
is a function defining the identifiability constraints of
the model.
Most Stochastic mortality models proposed in the literature can be cast to this representation (See Hunt and Blake (2015)).
Parametric age functions should be scalar functions of the form
f <- function(x, ages) taking a scalar age x and a vector
of model fitting ages (see examples below).
Do to limitation of functions gnm within package
gnm, which is used for fitting "StMoMo" objects to data
(see fit.StMoMo), models combining parametric and
non-parametric age-modulating functions are not supported at the moment.
A list with class "StMoMo" with components:
link |
a character string defining the link function of the model. |
staticAgeFun |
a logical value indicating if the model has a static age function. |
periodAgeFun |
a list defining the period age modulating parameters. |
cohortAgeFun |
an object defining the cohort age modulating parameters. |
constFun |
a function defining the identifiability constraints. |
N |
an integer specifying The number of age-period terms in the model. |
textFormula |
a character string of the model formula. |
gnmFormula |
a formula that can be used for fitting the model with package gnm. |
Plat, R. (2009). On stochastic mortality modeling. Insurance: Mathematics and Economics, 45(3), 393-404.
Hunt, A., & Blake, D. (2015). On the Structure and Classification of Mortality Models Mortality Models. Pension Institute Working Paper. http://www.pensions-institute.org/workingpapers/wp1506.pdf.
fit.StMoMo, lc, cbd,
apc, rh, m6, m7,
m8
#Lee-Carter model 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 = "log", staticAgeFun = TRUE, periodAgeFun = "NP", constFun = constLC) plot(fit(LC, data = EWMaleData, ages.fit = 55:89)) #CBD model f2 <- function(x, ages) x - mean(ages) CBD <- StMoMo(link = "logit", staticAgeFun = FALSE, periodAgeFun = c("1", f2)) plot(fit(CBD, data = EWMaleData, ages.fit = 55:89)) #Reduced Plat model (Plat, 2009) 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) #nsum g(c)=0, nsum cg(c)=0, nsum 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 #nsum 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 = "log", staticAgeFun = TRUE, periodAgeFun = c("1", f2), cohortAgeFun = "1", constFun = constPlat) plot(fit(PLAT, data = EWMaleData, ages.fit = 55:89)) #Models not supported ## Not run: MnotSup1 <- StMoMo(periodAgeFun = c(f2, "NP")) MnotSup1 <- StMoMo(periodAgeFun = f2, cohortAgeFun = "NP") ## End(Not run)#Lee-Carter model 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 = "log", staticAgeFun = TRUE, periodAgeFun = "NP", constFun = constLC) plot(fit(LC, data = EWMaleData, ages.fit = 55:89)) #CBD model f2 <- function(x, ages) x - mean(ages) CBD <- StMoMo(link = "logit", staticAgeFun = FALSE, periodAgeFun = c("1", f2)) plot(fit(CBD, data = EWMaleData, ages.fit = 55:89)) #Reduced Plat model (Plat, 2009) 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) #nsum g(c)=0, nsum cg(c)=0, nsum 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 #nsum 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 = "log", staticAgeFun = TRUE, periodAgeFun = c("1", f2), cohortAgeFun = "1", constFun = constPlat) plot(fit(PLAT, data = EWMaleData, ages.fit = 55:89)) #Models not supported ## Not run: MnotSup1 <- StMoMo(periodAgeFun = c(f2, "NP")) MnotSup1 <- StMoMo(periodAgeFun = f2, cohortAgeFun = "NP") ## End(Not run)
Create StMoMoData object suitable for fitting a Stochastic Mortality
Model using function fit.StMoMo.
StMoMoData(data, series = names(data$rate)[1], type = c("central", "initial"))StMoMoData(data, series = names(data$rate)[1], type = c("central", "initial"))
data |
demogdata object of type "mortality". It is either the
output from functions |
series |
name of series within |
type |
the type of exposure that should be included in the
output. The alternatives are |
A list with class "StMoMoData" with components:
Dxt |
matrix of deaths data. |
Ext |
matrix of observed exposures. |
ages |
vector of ages corresponding to rows of |
years |
vector of years corresponding to rows of |
type |
the type of exposure in the data. |
series |
name of the extracted series. |
label |
label of the data. |
## Not run: library(demography) NZdata <- hmd.mx(country = "NZL_NP", username = username, password = password, label = "New Zealand") NZStMoMo <- StMoMoData(NZdata, series = "male") summary(NZStMoMo) ## End(Not run)## Not run: library(demography) NZdata <- hmd.mx(country = "NZL_NP", username = username, password = password, label = "New Zealand") NZStMoMo <- StMoMoData(NZdata, series = "male") summary(NZStMoMo) ## End(Not run)