Constrained Covariance


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);