# Most efficient way to re-use correlation/covariance matrix multiple times

Hello!

I have a financial time series problem where I assume that the covariance matrix at any point is given by a combination of the asset volatilities and the correlation matrix, ie

Sigma[t] = quad_form_diag(correl_matrix, vol[t])


where correl_matrix is an [n.n] correlation matrix and vol is a length n vector.

I later have a likelihood at each point that is multivariate Student-t.

The likelihood function involves the term

(r[t] - mu)' * inverse(Sigma) * (r[t]-mu)
= (r[t] - mu)' * inverse(quad_form_diag(correl_matrix, vol[t])) * (r[t]-mu)


Obviously I don’t want to re-calculate inv(Sigma) (either implicitly or explicitly) at each time step due to the O(n^3) complexity.

Is there a recommended to repeatedly calculate likelihoods with these quadratic forms involved?

The algebraic way would be keeping the inverse of the correlation matrix and then re-using it repeatedly, eg

(r[t] - mu)' * quad_form_diag(inverse(correl_matrix), inverse(vol[t])) * (r[t]-mu)


…but I can appreciate that using the direct inverse might not be the most stable or efficient way to do it.

I know about the mdivide functions and could work with the lower triangular decomposition of correl_matrix, but I’m not sure whether repeatedly calling eg mdivide_left_tri_low(correl_matrix_L * vol[t], r[t] - mu) will incur a much higher operations count

Thanks!

3 Likes

I’m afraid there’s not an efficient way to do this simply with Stan built-ins. Intrinsically, there’s no way to get around the solve, as you need it for the determinant in the normalization term.

At least with Cholesky factors it’s more efficient and stable.

You might be able to push the volatility term out into a scaling term on the data and then post-adjust the normalizing factor.

\log \textrm{normal}(y | \mu, \sigma) = -log(\sigma) - \frac{1}{2} ((y - \mu) / \sigma)^2 + \textrm{const.}

An alternative way to do this is scale y and \mu by \sigma, as in

\log \textrm{normal}(y / \sigma | \mu / \sigma, 1) = -\frac{1}{2} ((y - \mu) / \sigma)^2 + \textrm{const.}

It’s only off by the log adjustment. I’m not the world’s best linear algebraist, but there might be some way to do something similar for the multi-t given that the determinant of a product is a product of the determinants and the determinant of a diagonal scaling term is trivial.

If not, then just rewriting the whole thing in terms of Cholesky factors should help.

1 Like