# Can one compute the autocorrelation of covariance matrices sampled by MCMC?

Imagine that we sample a covariance matrix from a Wishart distribution by MCMC.

At every iteration, we get a new sample matrix $S_i$ from the Wishart distribution.

Question: Given the trace that contains all the samples $S_1,...S_n$, can I plot the autocorrelation of these samples?

I have seen someone using the autocorrelation of $\log(\det(S))$, but I found no justification.

The way I see it, when drawing $p \times p$ matrices $S_1, S_2, \dots S_n$ from a Wishart distribution, you are really drawing $p\cdot (p+1)/2$ specifically related univariate random variables at each time step. I.e., only the $\text{vech}(S_i)$ part of any $S_i$ [ = the upper triangular part] is random, and symmetry gives you the rest. In other words, the autocorrelation is defined pairwise for any two entries of the $p\cdot (p+1)/2$-dimensional vectors $\text{vech}(S_i)$ and $\text{vech}(S_{i-h})$ for any $h>0$. Of course, this leaves you with a possibly infeasibly large number of $(p\cdot (p+1)/2)^2$ univariate autocorrelations to track, and since the entries of each $\text{vech}(S_i)$ will be meanfingully related to one another (see for example the definition given via draws from a normal here: https://en.wikipedia.org/wiki/Wishart_distribution), I could well imagine that you discard information by doing this univariate analysis. That being said, the univariate autocorrelations can be calculated entry-wise by first defining \begin{align} \bar{S} &= \frac{1}{n}\sum_{i=1}^n \text{vech}(S_i) \\ \bar{S_h} & = \frac{1}{n-h}\sum_{i=h+1}^n \text{vech}(S_i)\text{vech}(S_{i-h})^T. \end{align} Clearly, $\bar{S}$ is a natural estimator for the upper triangular part of the expectation (that you can replace with the true expectation of your Wishart distribution if it is known to you). Similarly, $\bar{S}_h$ is a natural estimator for the moment $\mathbb{E}(\text{vech}(S_i)\text{vech}(S_{i-h})^T)$. Lastly, noting that \begin{align} \text{Cov}(\text{vech}(S_i)\text{vech}(S_{i-h})^T) = \mathbb{E}(\text{vech}(S_i)\text{vech}(S_{i-h})^T) - \mathbb{E}(\text{vech}(S_i))\mathbb{E}(\text{vech}(S_i))^T, \end{align} one arrives at the autocorrelation estimates $A(h)$ for $\text{vech}(S_i)$ via \begin{align} A(h) &= \bar{S}_h - \bar{S}\bar{S}^T. \end{align} As mentioned before, this gives you the $(p\cdot (p+1)/2)^2$ autocorrelations of each Wishart matrix entry with each other Wishart matrix entry. If that is too much information to display, I think one strategy that you could take would be to define the univariate time series \begin{align} a(h) &= \frac{1}{(p\cdot (p+1)/2)^2}\sum_{i=1}^{(p\cdot (p+1)/2}\sum_{j=1}^{(p\cdot (p+1)/2}|A(h)_{ij}|, \end{align} i.e. you simply take the average of the absolute value of the autocorrelation. If you only care about positive autocorrelations and don't think that negative autocorrelation is detrimental, then spare yourself the absolute value. Similarly, if you think that autocorrelation along the diagonal is worse than off the diagonal or the other way around, you could add weights $w_{ij}$ that take this 'loss function' into account.