Hi!
Is there an efficient way of modeling the following VAR(1) model
\eta_{t} = B \eta_{t-1} + \epsilon_{t}, where \eta_{t} \in \mathbb{R}^{2} is assumed to be stationary. However, each component is defined with a unit variance, and the correlation between the outputs is of interest.
data {
int N_Y;
int N_T;
int N_eta;
matrix[N_T, N_Y] Y;
int<lower = 1, upper = N_T> time[N_T];
}
parameters {
vector[N_eta] eta_1;
vector[N_eta] eta1_Z[N];
real<lower = -1, upper = 1> beta_lag_mu;
matrix<lower = -1, upper = 1>[N_eta, N_eta] beta_lag;
cholesky_factor_corr[N_eta] L_phiY1;
}
transformed parameters {
corr_matrix[N_eta] eta_corr1 = multiply_lower_tri_self_transpose(L_phiY1);
matrix[N_eta, N_eta] zeta_cov = eta_corr1 - quad_form_sym(eta_corr1, beta_lag');
matrix[N_eta, N_eta] zeta_Lcov = cholesky_decompose(zeta_cov);
}
model {
vector[N_eta] mu_eta1[N]; //expected outcome for LVs
vector[N_eta] eta1[N];
for(n in 1:N){
eta1_Z[n] ~ std_normal();
if(time[n] == 1){
mu_eta1[n] = eta_1;
}
else if(time[n] > 1){
mu_eta1[n] = beta_lag*eta1[n-1];
}
eta1[n] = mu_eta1[n] + zeta_Lcov*eta1_Z[n];
}
Basically, I am assuming that the var(\epsilon_t) = \Sigma - B\Sigma B' with \Sigma as the correlation matrix of \eta_t. By using this approach, I’m getting a lot of
Exception: cholesky_decompose: Matrix m is not positive definite
even during active sampling.