I have a VAR model defined as:
Y_{1,it} = \alpha_{1} + \beta_{1i}*Y_{1,i(t-1)} + \beta_{2i}*Y_{2,i(t-1)} +\epsilon_{1,it}
Y_{2,it} = \alpha_{2} + \beta_{3i}*Y_{1,i(t-1)} + \beta_{4i}*Y_{2,i(t-1)} +\epsilon_{2,it}.
\beta_{ji} is defined as \beta_{ji} = \beta_{j} + U_{\beta ji}, for j = 1,2,3,4 where the random effects are defined as
\epsilon_{it} \sim MVN(0, \Sigma)
U_{\beta j} \sim MVN(0, \Phi).
I will also like to point out that \epsilon_{1,it} and \epsilon_{2,it} are correlated. Similarly, all the level-2 random effects U_{\beta j} are correlated.
In order for the model to be stationary, I simulated U_{\beta j} from a multivariate normal distribution with very small variances and no covariances. Therefore, my covariance matrix \Phi is a diagonal matrix, and my correlation matrix is almost an identity matrix.
Everything seems to work but I am having trouble recovering the parameter values for the correlation matrix for the level-2 random effects. I am not sure if I should use lkj_corr_cholesky(10) as prior for the random effects correlation matrix since some of my values are way off
My stan code is
data {
int N_Y;
int N_ID;
int N_T;
vector[2] Y[N_ID, N_T];
}
parameters {
vector[4] UB_I[N_ID]; //random effects
vector<lower = 0>[4] B_tau; //random effects SD
cholesky_factor_corr[4] B_ch; //cholesky correlated random effects
vector[4] Beta; //overall auto and cross lag regression coefficients
vector[2] alpha; //fixed intercepts
vector<lower = 0>[N_Y] Y_sd; //Level 1 residual sigma
cholesky_factor_corr[N_Y] Y_ch; //Level 1 residual cholesky corr
}
model {
vector[N_Y] mu_Y[N_ID,N_T]; //Expected Y
vector[4] B_I[N_ID]; //Level 2 auto and crossed lag slopes
for(j in 1: 4){
for(i in 1:N_ID){
B_I[i,j] = Beta[j] + UB_I[i,j];
}
}
for (i in 1:N_ID){
//time = 1
mu_Y[i,1,1] = alpha[1];
mu_Y[i,1,2] = alpha[2];
//time > 2
for(t in 2:N_T){
mu_Y[i,t,1] = alpha[1] + B_I[i,1]*Y[i,t-1,1] + B_I[i,2]*Y[i,t-1,2];
mu_Y[i,t,2] = alpha[2] + B_I[i,3]*Y[i,t-1,1] + B_I[i,4]*Y[i,t-1,2];
}
UB_I[i,] ~ multi_normal_cholesky(rep_vector(0,4), diag_pre_multiply(B_tau, B_ch));
for (t in 1:N_T){
Y[i,t] ~ multi_normal_cholesky(mu_Y[i,t], diag_pre_multiply(Y_sd, Y_ch));
}
}
//priors
Y_sd ~ normal(0,3);
Y_ch ~ lkj_corr_cholesky(1);
B_ch ~ lkj_corr_cholesky(10);
B_tau ~ std_normal();
Beta ~ std_normal();
}
generated quantities {
corr_matrix[4] B_corr;
cov_matrix[4] B_cov;
cov_matrix[N_Y] Y_cov;
corr_matrix[N_Y] Y_corr;
B_corr = multiply_lower_tri_self_transpose(B_ch);
B_cov = quad_form_diag(B_corr, B_tau);
Y_corr = multiply_lower_tri_self_transpose(Y_ch);
Y_cov = quad_form_diag(Y_corr, Y_sd);
}
My code for simulating the data is
N <- 200 #number of individuals
Nt <- 20 #Number of time occasions
eps <- vector("list", length = Nt) #Level-1 residual structure
Y1 <- vector("list", length = Nt)
#Residuals for Structural Model
sigma <- matrix(c(1, .2,.2,1),2)
for (j in 1:Nt) {
eps[[j]] <- mvtnorm::rmvnorm(N, mean = c(0,0), sigma = sigma)
}
#time 1
Y1 <- eps
cov_beta <- matrix(c(.01, 0, 0, 0,
0, .01, 0, 0,
0, 0, .02, 0,
0, 0, 0, .05), nrow = 4)
ov_beta <- c(.4, .2, .4,.2) #overall beta_lag
#autoregression slopes
beta_lag <- mvtnorm::rmvnorm(N, mean = ov_beta, sigma = cov_beta)
#time 2 & beyond
for (j in 2:Nt) {
Y1[[j]][,1] <- Y1[[j-1]][,1]*beta_lag[,1] + Y1[[j-1]][,2]*beta_lag[,2] + eps[[j]][,1]
Y1[[j]][,2] <- Y1[[j-1]][,1]*beta_lag[,3] + Y1[[j-1]][,2]*beta_lag[,4] + eps[[j]][,2]
}
Y <- array(0, dim = c(N,Nt,2))
for (j in 1:2) {
for (t in 1:Nt) {
for(i in 1:N) {
Y[i,t,j] = Y1[[t]][i,j]
}
}
}
