ex_brm.stan

// generated with brms 2.10.3
functions {
}
data {
  int<lower=1> N;  // number of observations
  vector[N] Y;  // response variable
  // data needed for ARMA correlations
  int<lower=0> Kar;  // AR order
  int<lower=0> Kma;  // MA order
  // number of lags per observation
  int<lower=0> J_lag[N];
  // data for group-level effects of ID 1
  int<lower=1> N_1;  // number of grouping levels
  int<lower=1> M_1;  // number of coefficients per level
  int<lower=1> J_1[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_1_1;
  vector[N] Z_1_2;
  vector[N] Z_1_3;
  vector[N] Z_1_sigma_4;
  vector[N] Z_1_sigma_5;
  int<lower=1> NC_1;  // number of group-level correlations
  // data for group-level effects of ID 2
  int<lower=1> N_2;  // number of grouping levels
  int<lower=1> M_2;  // number of coefficients per level
  int<lower=1> J_2[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_2_1;
  vector[N] Z_2_2;
  vector[N] Z_2_3;
  vector[N] Z_2_sigma_4;
  vector[N] Z_2_sigma_5;
  int<lower=1> NC_2;  // number of group-level correlations
  // data for group-level effects of ID 3
  int<lower=1> N_3;  // number of grouping levels
  int<lower=1> M_3;  // number of coefficients per level
  int<lower=1> J_3[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_3_1;
  vector[N] Z_3_2;
  vector[N] Z_3_3;
  vector[N] Z_3_sigma_4;
  vector[N] Z_3_sigma_5;
  int<lower=1> NC_3;  // number of group-level correlations
  int prior_only;  // should the likelihood be ignored?
}
transformed data {
  int max_lag = max(Kar, Kma);
  matrix[N, 3] Z_1;
  matrix[N, 2] Z_1_sigma;
  matrix[N, 3] Z_2;
  matrix[N, 2] Z_2_sigma;
  matrix[N, 3] Z_3;
  matrix[N, 2] Z_3_sigma;
  Z_1[, 1] = Z_1_1;
  Z_1[, 2] = Z_1_2;
  Z_1[, 3] = Z_1_3;

  Z_1_sigma[, 1] = Z_1_sigma_4;
  Z_1_sigma[, 2] = Z_1_sigma_5;

  Z_2[, 1] = Z_2_1;
  Z_2[, 2] = Z_2_2;
  Z_2[, 3] = Z_2_3;

  Z_2_sigma[, 1] = Z_2_sigma_4;
  Z_2_sigma[, 2] = Z_2_sigma_5;

  Z_3[, 1] = Z_3_1;
  Z_3[, 2] = Z_3_2;
  Z_3[, 3] = Z_3_3;

  Z_3_sigma[, 1] = Z_3_sigma_4;
  Z_3_sigma[, 2] = Z_3_sigma_5;

}
parameters {
  // temporary intercept for centered predictors
  real Intercept;
  // temporary intercept for centered predictors
  real Intercept_sigma;
  vector[Kar] ar;  // autoregressive effects
  vector[Kma] ma;  // moving-average effects
  vector<lower=0>[M_1] sd_1;  // group-level standard deviations
  matrix[M_1, N_1] z_1;  // standardized group-level effects
  // cholesky factor of correlation matrix
  cholesky_factor_corr[M_1] L_1;
  vector<lower=0>[M_2] sd_2;  // group-level standard deviations
  matrix[M_2, N_2] z_2;  // standardized group-level effects
  // cholesky factor of correlation matrix
  cholesky_factor_corr[M_2] L_2;
  vector<lower=0>[M_3] sd_3;  // group-level standard deviations
  matrix[M_3, N_3] z_3;  // standardized group-level effects
  // cholesky factor of correlation matrix
  cholesky_factor_corr[M_3] L_3;
}
transformed parameters {
  // actual group-level effects
  matrix[N_1, 3] r_1_mu;
  matrix[N_1, 2] r_1_sigma;
  matrix[N_2, 3] r_2_mu;
  matrix[N_2, 2] r_2_sigma;
  matrix[N_3, 3] r_3_mu;
  matrix[N_3, 2] r_3_sigma;

  matrix[N_1, M_1] r_1 = (diag_pre_multiply(sd_1, L_1) * z_1)';
  // actual group-level effects
  matrix[N_2, M_2] r_2 = (diag_pre_multiply(sd_2, L_2) * z_2)';
  // actual group-level effects
  matrix[N_3, M_3] r_3 = (diag_pre_multiply(sd_3, L_3) * z_3)';
  r_1_mu[, 1] = r_1[, 1];
  r_1_mu[, 2] = r_1[, 2];
  r_1_mu[, 3] = r_1[, 3];
  r_1_sigma[, 1] = r_1[, 4];
  r_1_sigma[, 2] = r_1[, 5];

  r_2_mu[, 1] = r_2[, 1];
  r_2_mu[, 2] = r_2[, 2];
  r_2_mu[, 3] = r_2[, 3];
  r_2_sigma[, 1] = r_2[, 4];
  r_2_sigma[, 2] = r_2[, 5];

  r_3_mu[, 1] = r_3[, 1];
  r_3_mu[, 2] = r_3[, 2];
  r_3_mu[, 3] = r_3[, 3];
  r_3_sigma[, 1] = r_3[, 4];
  r_3_sigma[, 2] = r_3[, 5];
}
model {
  // initialize linear predictor term
  vector[N] mu = Intercept + rep_vector(0, N);
  // initialize linear predictor term
  vector[N] sigma = Intercept_sigma + rep_vector(0, N);
  // objects storing residuals
  matrix[N, max_lag] Err = rep_matrix(0, N, max_lag);
  vector[N] err;
  mu += transpose(columns_dot_product(transpose(r_1_mu[J_1, ]), transpose(Z_1)));
  mu += transpose(columns_dot_product(transpose(r_2_mu[J_2, ]), transpose(Z_2)));
  mu += transpose(columns_dot_product(transpose(r_3_mu[J_3, ]), transpose(Z_3)));
  sigma += transpose(columns_dot_product(transpose(r_1_sigma[J_1, ]), transpose(Z_1_sigma)));
  sigma += transpose(columns_dot_product(transpose(r_2_sigma[J_2, ]), transpose(Z_2_sigma)));
  sigma += transpose(columns_dot_product(transpose(r_3_sigma[J_3, ]), transpose(Z_3_sigma)));
  // apply the inverse link function
  sigma = exp(sigma);
  // include ARMA terms
  for (n in 1:N) {
    mu[n] += Err[n, 1:Kma] * ma;
    err[n] = Y[n] - mu[n];
    for (i in 1:J_lag[n]) {
      Err[n + 1, i] = err[n + 1 - i];
    }
    mu[n] += Err[n, 1:Kar] * ar;
  }
  // priors including all constants
  target += normal_lpdf(Intercept | 0, 1);
  target += normal_lpdf(Intercept_sigma | 0, 1);
  target += normal_lpdf(ar | 0, 1.3)
    - 1 * log_diff_exp(normal_lcdf(1 | 0, 1.3), normal_lcdf(-1 | 0, 1.3));
  target += normal_lpdf(ma | 0, 1.3)
    - 1 * log_diff_exp(normal_lcdf(1 | 0, 1.3), normal_lcdf(-1 | 0, 1.3));
  target += normal_lpdf(sd_1[1] | 0.5, 1)
    - 1 * normal_lccdf(0 | 0.5, 1);
  target += student_t_lpdf(sd_1[2] | 5, 0.5, 1)
    - 1 * student_t_lccdf(0 | 5, 0.5, 1);
  target += student_t_lpdf(sd_1[3] | 5, 0.5, 1)
    - 1 * student_t_lccdf(0 | 5, 0.5, 1);
  target += normal_lpdf(sd_1[4] | 0.8, 1)
    - 1 * normal_lccdf(0 | 0.8, 1);
  target += student_t_lpdf(sd_1[5] | 5, 0.5, 1)
    - 1 * student_t_lccdf(0 | 5, 0.5, 1);
  target += normal_lpdf(to_vector(z_1) | 0, 1);
  target += lkj_corr_cholesky_lpdf(L_1 | 12);
  target += student_t_lpdf(sd_2[1] | 5, 1, 5)
    - 1 * student_t_lccdf(0 | 5, 1, 5);
  target += student_t_lpdf(sd_2[2] | 5, 0.5, 1)
    - 1 * student_t_lccdf(0 | 5, 0.5, 1);
  target += student_t_lpdf(sd_2[3] | 5, 0.5, 1)
    - 1 * student_t_lccdf(0 | 5, 0.5, 1);
  target += normal_lpdf(sd_2[4] | 0.5, 1)
    - 1 * normal_lccdf(0 | 0.5, 1);
  target += student_t_lpdf(sd_2[5] | 5, 0.5, 1)
    - 1 * student_t_lccdf(0 | 5, 0.5, 1);
  target += normal_lpdf(to_vector(z_2) | 0, 1);
  target += lkj_corr_cholesky_lpdf(L_2 | 12);
  target += student_t_lpdf(sd_3[1] | 5, 1, 5)
    - 1 * student_t_lccdf(0 | 5, 1, 5);
  target += student_t_lpdf(sd_3[2] | 5, 2, 5)
    - 1 * student_t_lccdf(0 | 5, 2, 5);
  target += student_t_lpdf(sd_3[3] | 5, 2, 5)
    - 1 * student_t_lccdf(0 | 5, 2, 5);
  target += normal_lpdf(sd_3[4] | 0.5, 1)
    - 1 * normal_lccdf(0 | 0.5, 1);
  target += student_t_lpdf(sd_3[5] | 5, 2, 5)
    - 1 * student_t_lccdf(0 | 5, 2, 5);
  target += normal_lpdf(to_vector(z_3) | 0, 1);
  target += lkj_corr_cholesky_lpdf(L_3 | 12);
  // likelihood including all constants
  if (!prior_only) {
    target += normal_lpdf(Y | mu, sigma);
  }
}
generated quantities {
  // actual population-level intercept
  real b_Intercept = Intercept;
  // actual population-level intercept
  real b_sigma_Intercept = Intercept_sigma;
  // group-level correlations
  corr_matrix[M_1] Cor_1 = multiply_lower_tri_self_transpose(L_1);
  vector<lower=-1,upper=1>[NC_1] cor_1;
  // group-level correlations
  corr_matrix[M_2] Cor_2 = multiply_lower_tri_self_transpose(L_2);
  vector<lower=-1,upper=1>[NC_2] cor_2;
  // group-level correlations
  corr_matrix[M_3] Cor_3 = multiply_lower_tri_self_transpose(L_3);
  vector<lower=-1,upper=1>[NC_3] cor_3;
  // extract upper diagonal of correlation matrix
  for (k in 1:M_1) {
    for (j in 1:(k - 1)) {
      cor_1[choose(k - 1, 2) + j] = Cor_1[j, k];
    }
  }
  // extract upper diagonal of correlation matrix
  for (k in 1:M_2) {
    for (j in 1:(k - 1)) {
      cor_2[choose(k - 1, 2) + j] = Cor_2[j, k];
    }
  }
  // extract upper diagonal of correlation matrix
  for (k in 1:M_3) {
    for (j in 1:(k - 1)) {
      cor_3[choose(k - 1, 2) + j] = Cor_3[j, k];
    }
  }
}