Covariance matrix error because of float precision

Modeling
1

Hi all.

I am trying to implement the BEKK model(https://arxiv.org/pdf/math/0702815.pdf) using Stan and I face the following problem. When adding terms like A * A^T it returns some exceptions (probably because of float precision). I would appreciate any help to avoid this problem. The following is that part of my code and the exception message.

Best,
Edgar

parameters {
    matrix[N,N] A;
    matrix[N,N] B;
    cov_matrix[N] S0;
}
transformed parameters {
    cov_matrix[N] S[T];t
    for(t in 1:T)
    {
        if (t == 1){
            S[t] = S0;
        }
        else{
            S[t] = S0 + A * y[t-1] * y[t-1]' * A' + B * S[t-1] * B';  // the line that 
                                                                      // throws the
                                                                      // exception
        }
    }
}

image
image.png1773×937 163 KB

2

You would to use functions like crossprod, tcrossprod, and quad_form to enforce symmetry. Something like

S[t] = S0 + tcrossprod(A * y[t-1]) + quad_form(S[t - 1], B');
1 Like
3

I’ve generally had to resort to passing my matrices through something like:

  matrix makesym(matrix mat, int verbose, int pd){
    matrix[rows(mat),cols(mat)] out;
    for(coli in 1:cols(mat)){
      if(pd ==1 && mat[coli,coli] < 1e-5){
        if(verbose > 0) print("diagonal too low (",mat[coli,coli],") during makesym row ", coli, " col ", coli);
        out[coli,coli] = 1e-5;
      } else out[coli,coli] = mat[coli,coli]; 
      for(rowi in coli:rows(mat)){
        if(rowi > coli) {
          out[rowi,coli] = mat[rowi,coli]; 
          out[coli,rowi] = mat[rowi,coli];
        }
      }
    }
    return out;
  }
1 Like