Hello all!

Im trying to do parameter estimation in a factor stochastic volatility model and I’m getting some strange behaviour when sampling. As in this question: Strange convergence behaviour, I can sometimes run the model and get good results and sometimes I get things like the plot below.

Any idea what may be the reason to this behaviour? Is is that, when it works, it does not explore the posterior good enough? In any case, is there anything one can do to get around this?

The code for the model is given below:

```
//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_std[t,]*sigma_x[t];
x[t,] = x[t,] + mu_x[t];
x[t,1] = x[t,1]/sqrt(1 - square(phi_x[t]));
//print("x[,t]:",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 = 10;
real a0 = 10;
real b0 = 3;
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));
mu_x ~ normal(b_mu,B_mu);
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);
}
}
generated quantities{
matrix[2,n] y_new;
matrix[2,2] Sigma; //Covariance matrix for observations as function of parameters
vector[2] null;
null[1] = 0;
null[2] = 0;
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_new[,i] = multi_normal_rng(null,Sigma);
}
}
```

I hope this makes sense. Thanks in advance!

Jens