# Multivariate factor stochastic volatility model - Problem with covariance matrix

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.

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
real<lower=0,upper=1> phi_h_std;
real<lower=0> sigma_h_std;
vector<lower=0,upper=1> phi_x_std;
vector mu_x;
vector<lower = 0> 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> phi_x;
vector<lower=0> 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 = h/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 null;
matrix[2,2] Sigma; //Covariance matrix for observations as function of parameters
null = 0;
null = 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)*exp(h[i]) + exp(x[1,i]);
Sigma[1,2] = beta*beta*exp(h[i]);
Sigma[2,2] = square(beta)*exp(h[i]) + exp(x[2,i]);
Sigma[2,1] = Sigma[1,2];
y[,i] ~ multi_normal(null,Sigma);
}
}

``````

If you haven’t done so already, try setting `init_r` to a value that is less than its default of 2 in order to have less extreme initial values.

Thanks for the response @bgoodri. I have tried different values of init_r, but without success. It’s strange, because I have build the same likelihood with the package TMB (template model builder) and used the package tmbstan to sample from this model in stan, but when I now tried to do it directly, it fails. The only difference is the reparameterization of the AR(1) processes.

I can’t see the problem either. I think you are just going to have to start doing stuff like

``````if (is_nan(Sigma[1,1]) || is_nan(Sigma[1,2]) || is_nan(Sigma[2,2])) {
print("h[i] = ", h[i]);
pint("x[,i] = ", x[,i]);
print("beta = ", beta);
print("beta = ", beta);
}
``````

until you see what the problem is.

Thank you!
I think this will help a lot, because now I get the error:

``````h[i] = 0.139476
x[,i] = [nan,nan]
beta = 0.489538
beta = 0.5904
``````

which shows that there is something wrong with the x process.

I found the error. Instead of multiplying

``````   x[t,] = x_std[t,]*sigma_x[t];
``````

``````    x[t,] = x[t,]*sigma_x[t];