Predict method for Stochastic Mortality Models fits
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"), ...)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
\hat{η}_{xt} = o_{xt} + α_x + ∑_{i=1}^N β_x^{(i)}κ_t^{(i)} + β_x^{(0)}γ_{t-x}
for a fitted Stochastic Mortality model.
In producing a prediction the static age function, α_x, and the
age-modulating parameters, β_x^{(i)}, i=0, ..., N, are taken from
the fitted model in object while the period indexes,
κ_t^{(i)}, i=1,..., N, and cohort index, γ_{t-x},
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)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.