# Constrained Covariance

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.

Try adding a small epsilon to the diagonal. See for more info Cholesky decomposition failed for precision/covariance matrix - #5 by anon75146577

...
matrix[N_eta, N_eta] zeta_Lcov = add_diag(cholesky_decompose(zeta_cov), 1e-6);