Way back @bgoodri was referenced in @James_Savage post about hierarchical vars at https://khakieconomics.github.io/2016/11/27/Hierarchical-priors-for-panel-vector-autoregressions.html, concerning estimation of a weighted correlation matrix between global and local correlation matrices,
with
My question concerns using the cholesky corr formulation. So instead of \Omega I have it’s cholesky corr (ie I have L from \Omega = \text{LL}'). Now it obviously (or maybe not so obviously) is not equivalent to just weight the corresponding cholesky global and local correlation matrices with \rho. So it means I’m reconstructing the full correlation matrices, doing the weighting and then decomping again,
Omega[i] = cholesky_decompose(rho * multiply_lower_tri_self_transpose(Omega_global) + (1 - rho) * multiply_lower_tri_self_transpose(Omega_local[i]))
Is it possible to update the known cholesky factors and avoid the decomposition?
Some more info on the error when you just weight the cholesky corr matrices:
Chain 1 Exception: hier_model_namespace::log_prob: Omega[sym1__] is not a valid unit vector. The sum of the squares of the elements should be 1, but is 0.999903
and the stan model with cholesky corrs
data {
int N; // number of observations (number of rows in the panel)
int K; // number of dimensions of the outcome variable
int I; // number of individuals
int T; // The greatest number of time periods for any individual
int<lower = 1, upper = I> individual[N]; // integer vector the same length
// as the panel, indicating which individual each row corresponds to
int<lower = 1, upper = T> time[N]; // integer vector with individual-time
// (not absolute time). Each individual starts at 1
matrix[N, K] Y; // the outcome matrix, each individual stacked on each
// other, observations ordered earliest to latest
}
transformed data{
matrix[N - 1, K] Y_lag = Y[2:N];
}
parameters {
// individual-level parameters
cholesky_factor_corr[K] L_Omega_local[I]; // cholesky factor of correlation matrix of
// residuals each individual (across K outcomes)
vector<lower = 0>[K] tau[I]; // scale vector for residuals
matrix[K, K] z_beta[I];
vector[K] z_alpha[I];
// hierarchical priors (individual parameters distributed around these)
real<lower = 0, upper = 1> rho;
cholesky_factor_corr[K] L_Omega_global;
vector[K] tau_location;
vector<lower =0>[K] tau_scale;
matrix[K,K] beta_hat_location;
matrix<lower = 0>[K,K] beta_hat_scale;
vector[K] alpha_hat_location;
vector<lower = 0>[K] alpha_hat_scale;
}
transformed parameters {
matrix[K, K] beta[I]; // VAR(1) coefficient matrix
vector[K] alpha[I]; // intercept matrix
cholesky_factor_corr[K] L_Omega[I];
for(i in 1:I) {
alpha[i] = alpha_hat_location + alpha_hat_scale .*z_alpha[i];
beta[i] = beta_hat_location + beta_hat_scale .*z_beta[i];
L_Omega[i] = cholesky_decompose(rho * multiply_lower_tri_self_transpose(L_Omega_global) + (1 - rho) * multiply_lower_tri_self_transpose(L_Omega_local[i]));
}
}
model {
// hyperpriors
rho ~ beta(2, 2);
tau_location ~ cauchy(0, 1);
tau_scale ~ cauchy(0, 1);
alpha_hat_location ~ normal(0, 1);
alpha_hat_scale ~ cauchy(0, 1);
to_vector(beta_hat_location) ~ normal(0, .5);
to_vector(beta_hat_scale) ~ cauchy(0, .5);
L_Omega_global ~ lkj_corr_cholesky(1);
// hierarchical priors
for(i in 1:I) {
// non-centered parameterization
z_alpha[i] ~ normal(0, 1);
to_vector(z_beta[i]) ~ normal(0, 1);
tau[i] ~ normal(tau_location, tau_scale);
L_Omega_local[i] ~ lkj_corr_cholesky(10);
}
// likelihood
Y[2:N] ~ multi_normal_cholesky(alpha[individual] + beta[individual] * Y_lag',
diag_pre_multiply(tau[individual], L_Omega[individual]));
}