Hello,

I’m trying to fit a 2 dimensional factor stochastic volatility model:

y_1(t) = beta_1 * exp(h(t)) * epsilon(t) + exp(x_1(t))*eta_1(t)

y_2(t) = beta_2 * exp(h(t)) * epsilon(t) + exp(x_1(t)) * eta_2(t),

where epsilon, eta_1 and eta_2 are iid standard normal, h(t) is a centered AR(1) process and x_1(t) and x_2(t) are non-centered AR(1) processes, all independent.

I’m getting an error message concerning the covariance matrix:

```
Rejecting initial value:
Error evaluating the log probability at the initial value.
Exception: multi_normal_lpdf: LDLT_Factor of covariance parameter is not positive definite. last conditional variance is nan. (in 'model29277df03a02_2DimFactor' at line 85)
```

I have also tried to do a Cholesky decomposition of the covariance matrix and used multi_normal_cholesky, but without success.

Any tips on why this does not work? I have earlier tried to do this in a likelihood framework by using the Laplace approximation, but it is very unstable and I’m therefore trying Stan.

Thanks in advance!

Jens

```
//2 dim factor model
data {
int<lower=0> n; // Sample size
matrix[2,n] y; // data
}
parameters {
vector[n] h_std; //standardized AR(1) for factor
matrix[2,n] x_std; //standardized AR(1) idiosyncratic process
vector[2] beta; //loading matrix (vector in this case)
real<lower=0,upper=1> phi_h_std;
real<lower=0> sigma_h_std;
vector<lower=0,upper=1>[2] phi_x_std;
vector[2] mu_x;
vector<lower = 0>[2] sigma_x_std;
}
transformed parameters {
matrix[2,n] x; //indiosyncratic processes
vector[n] h; // factor volatility
real<lower=-1,upper=1> phi_h;
real<lower=0> sigma_h;
vector<lower=-1,upper=1>[2] phi_x;
vector<lower=0>[2] sigma_x;
phi_h = 2*phi_h_std - 1;
sigma_h = sqrt(sigma_h_std);
phi_x = 2*phi_x_std - 1;
sigma_x = sqrt(sigma_x_std);
//Scale with standard deviation
//Stationary assumption on h
h = h_std * sigma_h;
h[1] = h[1]/sqrt(1 - square(phi_h));
//Scale with standard deviation
//Stationary assumption on x
for(t in 1:2){
x[t,] = x[t,]*sigma_x[t];
x[t,] = x[t,] + mu_x[t];
x[t,1] = x[t,1]/sqrt(1 - square(phi_x[t]));
}
for(t in 2:n){
h[t] = h[t] + phi_h*h[t-1];
for(i in 1:2){
x[i,t] = x[i,t] + phi_x[i]*(x[i,t-1] - mu_x[i]);
}
}
}
model {
//Prior parameters
real B_beta = 1;
real b_mu = 0;
real B_mu = 100;
real a0 = 20;
real b0 = 1.5;
real B_sigma = 1;
vector[2] null;
matrix[2,2] Sigma; //Covariance matrix for observations as function of parameters
null[1] = 0;
null[2] = 0;
//Priors
sigma_h_std ~ gamma(0.5,1/(2*B_sigma));
sigma_x_std ~ gamma(0.5,1/(2*B_sigma));
phi_h_std ~ beta(a0,b0);
phi_x_std ~ beta(a0,b0);
beta ~ normal(null,B_beta);
h_std ~ normal(0,1);
x_std[1,] ~ normal(0,1);
x_std[2,] ~ normal(0,1);
//Constribution from observations
for(i in 1:n){
Sigma[1,1] = square(beta[1])*exp(h[i]) + exp(x[1,i]);
Sigma[1,2] = beta[1]*beta[2]*exp(h[i]);
Sigma[2,2] = square(beta[2])*exp(h[i]) + exp(x[2,i]);
Sigma[2,1] = Sigma[1,2];
y[,i] ~ multi_normal(null,Sigma);
}
}
```