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.

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

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[1] = ", beta[1]);
  print("beta[2] = ", beta[2]);
}

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[1] = 0.489538
beta[2] = 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];

I had

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

It was a easy to find when using the code you provided.

Thank you very much. I really appreciate the help!

Jens