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"), ... )
| 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. |
| ... | Arguments passed to |
An object of class stanfit returned by rstan::sampling.
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).$$
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.
#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: