How to test if a cross-covariance matrix is non-zero? The background of my study: 
In a Gibbs sampling where we sample $X$ (the variable of interests) and $Y$ from $P(X|Y)$ and $P(Y|X)$ respectively, where $X$ and $Y$ are $k$-dimensional random vectors. We know that the process is usually split into two stages:


*

*Burn-in Period, where we discard all the samples. Denote the samples as $X_1\sim X_t$ and $Y_1\sim Y_t$.

*"After-Burn-in" Period, where we average the samples $\bar{X} = \frac{1}{k}\sum_{i=1}^k X_{t+i}$ as our final desired result.


However, the samples in the "after-burn-in" sequence $X_{t+1}\sim X_{t+k}$ are not independently distributed. Therefore if I want to inspect the variance of the final result, it becomes
$$\operatorname{Var}[\bar{X}] = \operatorname{Var}\left[\sum_{i=1}^k X_{t+i}\right] = \frac{1}{k^2}\left(\sum_{i=1}^k\operatorname{Var}[X_{t+i}] + \sum_{i=1}^{k-1} \sum_{j=i+1}^k \operatorname{Cov}[X_{t+i},X_{t+j}]\right)$$
Here the term $\operatorname{Cov}[X_{t+i},X_{t+j}]$ is a $k\times k$ cross-covariance matrix applies to any $(i,j)$ with $i<j$. 
For example, I have 
$$X_{t+1} = (1,2,1)'\\
X_{t+2} = (1,0,2)'\\
X_{t+3} = (1,0,0)'\\
X_{t+4} = (5,0,-1)'$$
then I could estimate the covariance matrix $\operatorname{Cov}[X_{t+i}, X_{t+i+1}]$ with 
$$
\frac{1}{3}\sum_{i=1}^3 (X_{t+i}-\mu_{t+i})(X_{t+i+1}-\mu_{t+i+1})'
$$
Now I am interested in if the resulting estimation is significantly non-zero so that I need to include it into my variance estimation of $\operatorname{Var}[\bar{X}]$.
So here comes my questions:


*

*We sample $X_{t+i}$ from $P(X_{t+i}|Y_{t+i})$. Since $Y_{t+i}$ is changing, I think $X_{t+i}$ and $X_{t+i+1}$ are not from the same distribution, so $\operatorname{Cov}[X_{t+i},X_{t+j}]$ is not the same as $\operatorname{Cov}[X_{t+i},X_{t+i}]$. Is this statement correct?

*Suppose I have enough data to estimate $\operatorname{Cov}[X_{t+i},X_{t+i+1}]$ (neighboring samples in the sequence), is there any way to test if the covariance matrix is significantly a non-zero matrix? Broadly speaking, I am interested in an indicator which guides me to some meaningful cross-covariance matrices that should be included in my final variance estimation.

 A: 
  
*
  
*We sample $X_{t+i}$ from $P(X_{t+i}|Y_{t+i})$. Since $Y_{t+i}$ is changing, I think $X_{t+i}$ and $X_{t+i+1}$ are not from the same distribution [...]
  

You are confusing conditional and unconditional distributions here, see also my next remark. Conditional on $Y_{t+i} = y_1$ and $Y_{t+i+1} = y_2$, $P(X_{t+i}|Y_{t+i} = y_1) \neq P(X_{t+i+1}|Y_{t+i+1} = y_2)$. But the entire point of constructing your Gibbs sampler is to sample from the stationary distributions of $X$ and $Y$. Very roughly speaking, if you have run your chain for long enough and so that $\{Y_t\}$ follows the stationary distribution, you can then say
\begin{align}
 P(X_t) = \int_{\mathcal{Y}}P(X_t|Y_t)dP(Y_t),
\end{align}
meaning that the unconditional distribution of $X_t$ is also invariant. In other words, as $t \to \infty$ and we converge to the stationary distributions, $P(X_{t+i}|Y_{t+i}) = P(X_{t+i+1}|Y_{t+i+1})$, since $Y_{t+i}$ and $Y_{t+i+1}$ will asymptotically be drawn from (the same!) stationary distribution $P(Y_t)$. On the other hand and as before, once we condition on $Y_{t+i} = y_1$ and $Y_{t+i+1} = y_2$, this won't hold anymore, regardless how large $t$ is.

[...] so $\operatorname{Cov}[X_{t+i},X_{t+j}]$ is not the same as $\operatorname{Cov}[X_{t+i},X_{t+i}]$. Is this statement correct?

Yes, this is correct - even though $X_{t+1} \sim X_{t}$, i.e. $X_t$ and $X_{t+1}$ have the same stationary distribution. I know this may be confusing, but bear with me. Define $Y_t = 0.8\cdot Y_{t-1} + \varepsilon_t$ with $\varepsilon_t \overset{iid}{\sim} N(0,1)$. By iterated substitution, one can show that $Y_t = \sum_{i=0}^t0.8^i \varepsilon_{t-i}$, and since (infinite) sums of normals are still normal, it holds that $\text{Var}(Y_t) = \sum_{i=0}^t0.8^{2i} = \dfrac{1}{1-0.8^2}$ and so that $Y_t \overset{iid}{\sim} N(0, \dfrac{1}{1-0.8^2})$. Clearly, $Y_t$ and $Y_{t+1}$ will still be correlated, but they will also come from the same distribution ($Y_{t+1} \sim Y_{t}$). A similar situation holds for your $X_t$.


  
*Suppose I have enough data to estimate $\operatorname{Cov}[X_{t+i},X_{t+i+1}]$ (neighboring samples in the sequence), is there any way to test if the covariance matrix is significantly a non-zero matrix? Broadly speaking, I am interested in an indicator which guides me to some meaningful cross-covariance matrices that should be included in my final variance estimation.
  

Well, if you had infinitely many observations, they will all be significant eventually. Clearly, you cannot do this in practice, but there are ways of 'chopping off' the expansion after some terms, see the accepted excellent answer here. Basically, you define a kernel $k(\cdot)$ which decays to $0$ and assigns weights to the first $l_T$ covariance matrices that you could compute. If you want to choose $l_T$ in a principled way, you will have to dig a bit into the literature, but the post I linked gives you some good references to do exactly that.
