# R/simulatebootStMoMo.R In StMoMo: Stochastic Mortality Modelling

#### Documented in simulate.bootStMoMo

#' Simulate future sample paths from a Bootstrapped Stochastic
#' Mortality Model
#'
#' Simulate future sample paths from a Bootstrapped Stochastic Mortality Model.
#' The period indexes \eqn{\kappa_t^{(i)}, i = 1,..N,} are modelled
#' using ether a Multivariate Random Walk with Drift (MRWD) or
#' \eqn{N} independent ARIMA\eqn{(p, d, q)} models. The cohort index
#' \eqn{\gamma_{t-x}} is modelled using an ARIMA\eqn{(p, d, q)}.
#' By default an ARIMA\eqn{(1, 1, 0)} with a constant is used.
#'
#' @param object an object of class \code{"bootStMoMo"} with the bootstrapped
#' parameters of a stochastic mortality model.
#' @param nsim number of sample paths to simulate from each bootstrapped
#' sample. Thus if there are \code{nBoot} bootstrapped samples the total
#' number of paths will be \code{nsim * nBoot}.
#' @inheritParams simulate.mrwd
#' @param 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.
#' @inheritParams forecast.fitStMoMo
#'
#' @return A list of class \code{"simStMoMo"} with components
#'
#' \item{rates}{ a three dimensional array with the future simulated rates.}
#'
#' \item{ages}{ vector of ages corresponding to the first dimension of
#' \code{rates}.}
#'
#' \item{years}{ vector of years for which a simulations has been produced.
#' This corresponds to the second dimension of \code{rates}.}
#'
#' \item{kt.s}{ information on the simulated paths of the period indices of
#' the model. This is a list with the simulated paths of \eqn{\kappa_t}
#'  (\code{sim}) and the \code{years} for which simulations were produced.
#'  If the mortality model does not have any age-period terms (i.e. \eqn{N=0})
#'  this is set to \code{NULL}.}
#'
#' \item{gc.s}{ information on the simulated paths of the cohort index of the
#'  model. This is a list with the simulated paths of \eqn{\gamma_c}
#'  (\code{sim}) and the \code{cohorts} for which simulations were produced.
#'  If the mortality model does not have a cohort effect this is set to
#'  \code{NULL}.}
#'
#' \item{oxt.s}{ a three dimensional array with the offset used in the
#' simulations.}
#'
#' \item{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).}
#'
#'  \item{jumpchoice}{Jump-off method used in the simulation.}
#'
#' \item{kt.method}{method used in the modelling of the period index.}
#'
#'  \item{model}{the bootstrapped model from which the simulations were
#'  produced.}
#'
#' @details
#' For further details see \code{\link{simulate.fitStMoMo}}.
#'
#'
#' @examples
#' #Long computing times
#' \dontrun{
#' #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) #' } #' #' @export simulate.bootStMoMo <-function(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, ...) { jumpchoice <- match.arg(jumpchoice) kt.method <- match.arg(kt.method) ages <- object$model$ages nAges <- length(ages) nBoot <- length(object$bootParameters)
nPath <- nsim * nBoot

## Handle generator seed
if (!exists(".Random.seed", envir = .GlobalEnv))
runif(1)
if (is.null(seed))
RNGstate <- .Random.seed
else {
R.seed <- .Random.seed
set.seed(seed)
RNGstate <- structure(seed, kind = as.list(RNGkind()))
on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
}

## Offset
if (is.null(oxt))
oxt <- 0
oxt.s <- array(oxt, dim = c(nAges, h, nBoot))

## Do simulations

#Initialse output variables
modelFit <- object$bootParameters[[1]] modelFit$Dxt <- object$model$Dxt
modelFit$Ext <- object$model$Ext tempSim <- simulate.fitStMoMo(object = modelFit, nsim = nsim, seed = NULL, h = h, oxt = oxt.s[, , 1], gc.order = gc.order, gc.include.constant = gc.include.constant, jumpchoice = jumpchoice, kt.method = kt.method, kt.order = kt.order, kt.include.constant = kt.include.constant, kt.lookback = kt.lookback, gc.lookback = gc.lookback, ...) rates <- array(NA, c(dim(tempSim$rates)[1:2], nPath),
list(dimnames(tempSim$rates)[[1]], dimnames(tempSim$rates)[[2]], 1:nPath))
fitted <- array(NA, c(dim(tempSim$fitted)[1:2], nPath), list(dimnames(tempSim$fitted)[[1]],
dimnames(tempSim$fitted)[[2]], 1:nPath)) if (is.null(tempSim$kt.s)) {
kt.s <- NULL
} else {
kt.s <- list(sim = NULL, years = tempSim$kt.s$years)
kt.s$sim <- array(NA, c(dim(tempSim$kt.s$sim)[1:2], nPath), list(dimnames(tempSim$kt.s$sim)[[1]], dimnames(tempSim$kt.s$sim)[[2]], 1:nPath)) } if (is.null(tempSim$gc.s)) {
gc.s <- NULL
} else {
gc.s <- list(sim = NULL, cohorts = tempSim$gc.s$cohorts)
gc.s$sim <- array(NA, c(dim(tempSim$gc.s$sim)[1], nPath), list(dimnames(tempSim$gc.s$sim)[[1]], 1:nPath)) } # Do simulations for each bootstrap sample i <- 1 while (i <= nBoot) { rates[, , ((i - 1) * nsim + 1):(i * nsim)] <- tempSim$rates
fitted[, , ((i - 1) * nsim + 1):(i * nsim)] <- tempSim$fitted if (!is.null(kt.s)) kt.s$sim[, , ((i - 1) * nsim + 1):(i * nsim)] <- tempSim$kt.s$sim
if (!is.null(gc.s))
gc.s$sim[, ((i - 1) * nsim + 1):(i * nsim)] <- tempSim$gc.s$sim i <- i + 1 if ( i <= nBoot) { modelFit <- object$bootParameters[[i]]
modelFit$Dxt <- object$model$Dxt modelFit$Ext <- object$model$Ext
tempSim <- simulate.fitStMoMo(object = modelFit,
nsim = nsim, seed = NULL, h = h,
oxt = oxt.s[, , i], gc.order = gc.order,
gc.include.constant = gc.include.constant,
jumpchoice = jumpchoice,
kt.method = kt.method,
kt.order = kt.order,
kt.include.constant = kt.include.constant,
kt.lookback = kt.lookback,
gc.lookback = gc.lookback, ...)

}
}
structure(list(rates = rates, ages = tempSim$ages, years = tempSim$years, kt.s = kt.s, gc.s = gc.s,
oxt.s = oxt.s, fitted = fitted, jumpchoice = jumpchoice,
kt.method = kt.method, model = object, call = match.call()),
class = "simStMoMo")
}


## Try the StMoMo package in your browser

Any scripts or data that you put into this service are public.

StMoMo documentation built on May 2, 2019, 11:42 a.m.