Fit and Forecast Bayesian M6 model (CBD with cohort effect) introduced in Cairns et al (2009). The model can be fitted with a Poisson or Negative-Binomial distribution. The function outputs posteriors distributions for each parameter, predicted death rates and log-likelihoods.

m6_stan(
  death,
  exposure,
  forecast,
  age,
  validation = 0,
  family = c("poisson", "nb"),
  ...
)

Arguments

death

Matrix of deaths.

exposure

Matrix of exposures.

forecast

Number of years to forecast.

age

Vector of ages.

validation

Number of years for validation.

family

specifies the random component of the mortality model. "Poisson" assumes a Poisson model with log link and "nb" assumes a negative-binomial model with log link and overdispersion parameter \(\phi\).

...

Arguments passed to rstan::sampling (e.g. iter, chains).

Value

An object of class stanfit returned by rstan::sampling.

Details

The created model is either a log-Poisson or a log-Negative-Binomial version of the M6 model: $$D_{x,t} \sim \mathcal{P}(\mu_{x,t} e_{x,t})$$ or $$D_{x,t}\sim NB\left(\mu_{x,t} e_{x,t},\phi\right)$$ with $$\log \mu_{xt} = \kappa_t^{(1)} + (x-\bar{x})\kappa_t^{(2)}+\gamma_{t-x},$$ where \(\bar{x}\) is the average age in the data.

To ensure the identifiability of th model, we impose $$\gamma_1=0,\gamma_C=0,$$ where \(C\) represents the most recent cohort in the data.

For the period terms, we consider a multivariate random walk with drift: $$\boldsymbol{\kappa}_{t}=\boldsymbol{c}+ \boldsymbol{\kappa}_{t-1}+\boldsymbol{\epsilon}_{t}^{\kappa},\quad \boldsymbol{\kappa}_{t}=\left(\begin{array}{c}\kappa_{t}^{(1)} \\\kappa_{t}^{(2)}\end{array}\right), \quad \boldsymbol{\epsilon}_{t}^{\kappa} \sim N\left(\mathbf{0}, \Sigma\right),$$ with normal priors: \(\boldsymbol{c} \sim N(0,10)\).

The variance-covariance matrix of the error term is defined by $$\boldsymbol{\Sigma}=\left(\begin{array}{cc}\sigma_1^{2} & \rho_{\Sigma} \sigma_1 \sigma_2 \\\rho_{\Sigma} \sigma_1 \sigma_{Y} & \sigma_2^{2}\end{array}\right)$$ where the variance coefficients have independent exponential priors: \(\sigma_1, \sigma_2 \sim Exp(0.1)\) and the correlation parameter has a uniform prior: \(\rho_{\Sigma} \sim U\left[-1,1\right]\). As for the other models, the overdispersion parameter has a prior distribution given by $$\frac{1}{\phi} \sim Half-N(0,1).$$

For the cohort term, we consider a second order autoregressive process (AR(2)): $$\gamma_{c}=\psi_1 \gamma_{c-1}+\psi_2 \gamma_{c-2}+\epsilon^{\gamma}_{t},\quad \epsilon^{\gamma}_{t}\sim N(0,\sigma_{\gamma}).$$

To close the model specification, we impose some vague priors assumptions on the hyperparameters: $$\psi_1,\psi_2 \sim N(0,10),\quad \sigma_{\gamma}\sim Exp(0.1).$$

References

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.

Examples



#10-year forecasts for French data for ages 50-90 and years 1970-2017 with a log-Poisson model
ages.fit<-70:90
years.fit<-1990:2010
deathFR<-FRMaleData$Dxt[formatC(ages.fit),formatC(years.fit)]
exposureFR<-FRMaleData$Ext[formatC(ages.fit),formatC(years.fit)]
iterations<-50 # Toy example, consider at least 2000 iterations
fitM6=m6_stan(death = deathFR,exposure=exposureFR, age=ages.fit,forecast = 5,
family = "poisson",iter=iterations,chains=1)
#> 
#> SAMPLING FOR MODEL 'M6model' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.004 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 40 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: WARNING: There aren't enough warmup iterations to fit the
#> Chain 1:          three stages of adaptation as currently configured.
#> Chain 1:          Reducing each adaptation stage to 15%/75%/10% of
#> Chain 1:          the given number of warmup iterations:
#> Chain 1:            init_buffer = 3
#> Chain 1:            adapt_window = 20
#> Chain 1:            term_buffer = 2
#> Chain 1: 
#> Chain 1: Iteration:  1 / 50 [  2%]  (Warmup)
#> Chain 1: Iteration:  5 / 50 [ 10%]  (Warmup)
#> Chain 1: Iteration: 10 / 50 [ 20%]  (Warmup)
#> Chain 1: Iteration: 15 / 50 [ 30%]  (Warmup)
#> Chain 1: Iteration: 20 / 50 [ 40%]  (Warmup)
#> Chain 1: Iteration: 25 / 50 [ 50%]  (Warmup)
#> Chain 1: Iteration: 26 / 50 [ 52%]  (Sampling)
#> Chain 1: Iteration: 30 / 50 [ 60%]  (Sampling)
#> Chain 1: Iteration: 35 / 50 [ 70%]  (Sampling)
#> Chain 1: Iteration: 40 / 50 [ 80%]  (Sampling)
#> Chain 1: Iteration: 45 / 50 [ 90%]  (Sampling)
#> Chain 1: Iteration: 50 / 50 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.021 seconds (Warm-up)
#> Chain 1:                0.027 seconds (Sampling)
#> Chain 1:                0.048 seconds (Total)
#> Chain 1: