data { int J; // number of age categories int T; // number of years int d[J*T]; // vector of deaths vector[J* T] e; // vector of exposures vector[J] age; // vector of ages int family; // if family = 1 Neg-binomial } transformed data { real ma = mean(age); } parameters { real aux[family]; // neg. binomial dispersion parameter row_vector[2] c1; // parameters for kappa_1 row_vector[2] delta; // parameters for kappa_2 vector[2] rho; // autoregresive parameter real sigma[2]; // standard deviations for kappa_1 and kappa_2 real rho_sig; // correlation parameter matrix[T,2] epsilon; // *residual mean } transformed parameters { cov_matrix[2] Sig; // epsilon correlation matrix Sig[1,1] = sigma[1]; Sig[2,2] = sigma[2]; Sig[1,2] = prod(sigma)*rho_sig; } model { vector[J * T] mu; // process mean matrix[T,2] k; // process filter int pos = 1; // weird position parameter for (t in 1:T){ for (x in 1:J) { mu[pos] = log(e[pos]) + k[t,1]+(age[x]-ma)*k[t,2]; //Predictor dynamics pos += 1; } if(t == 1) k[1] = c1 + delta; else{ k[t] = c1 + delta*diag_matrix(rep_vector(t,2)); k[t] += k[t-1]*diag_matrix(rho) + epsilon[t]; } target += multi_normal_lpdf(epsilon[t]|rep_vector(0,2),Sig); } // Likelihood if (family == 0){ // Poisson log model target += poisson_log_lpmf(d |mu); } else { // Negative-Binomial log model target +=neg_binomial_2_log_lpmf (d|mu,inv(aux[1])); target += normal_lpdf(aux[1]|0,10); } // Priors target += exponential_lpdf(sigma | 0.1); target += normal_lpdf(c1|0,sqrt(10)); target += normal_lpdf(delta|0,sqrt(10)); target += normal_lpdf(rho|0,sqrt(10)); target += normal_lpdf(rho|0,sqrt(10)); }