Covariance of sum of two dependent random vectors

Question

$\newcommand{\cov}{\operatorname{cov}}\newcommand{\corr}{\operatorname{corr}}$I have two dependent multivariate random variables, or random vectors, $X$ and $Y$, their respective covariance matrices $\Sigma_X = \cov(X,X)$, $\Sigma_Y = \cov(Y,Y)$ and a non-zero correlation matrix $P_{XY} = \corr(X, Y)$.

I would like to compute the covariance matrix of the sum of the two random vectors:

\begin{align} Z {}={}& X + Y \\ \Sigma_Z {}={}& \Sigma_X + \Sigma_Y + 2 \cov(X, Y) \end{align}

Intuitively, I would assume the cross-covariance can be computed as follows:

$$\cov(X, Y) = P_{XY} \sqrt{\Sigma_X} \sqrt{\Sigma_Y}$$

But I can see that this matrix is not necessarily symmetric, and therefore the resulting $\Sigma_Z$ wouldn't be either. Am I wrong? Is it at all possible to compute $\Sigma_Z$ from the information I have? If not, is there some way to approximate it?

EDIT: I think it isn't correct to write:

$$ \Sigma_Z = \Sigma_X + \Sigma_Y + 2 \cov(X, Y) $$

as $\cov(X,Y) \ne \cov(Y,X)$, but instead $\cov(X,Y) = \cov(Y,X)^T$.

Rewriting the equation for $\Sigma_Z$ to:

\begin{align} \Sigma_Z {}={}& \Sigma_X + \Sigma_Y + \cov(X, Y) + \cov(Y, X) \\ {}={}& \Sigma_X + \Sigma_Y + \cov(X, Y) + \cov(X, Y)^T \end{align}

I can now prove that indeed $\Sigma_Z$ is guaranteed to be symmetric. As antkam points out in the chat below, it isn't necessary that $\cov(X, Y) + \cov(X, Y)^T$ be positive semi-definite for $\Sigma_Z$ to be positive semi-definite.

Answer 1

Your original intuition does not work, i.e., $cov(X, Y) \neq P_{XY} \sqrt{\Sigma_X} \sqrt{\Sigma_Y}$, because the definition of $corr(X_i, Y_j)$ involves scalar divisions, which cannot be reversed like this. For one thing, you have three matrices, so how do you pick the order of multiplication?

Instead, I would go back to the per-entry definitions:

$$[cov(X,Y)]_{ij} = cov(X_i, Y_j) = corr(X_i, Y_j) \ \sigma(X_i) \ \sigma(Y_j) = [P_{XY}]_{ij} \sqrt{ [\Sigma_X]_{ii}} \sqrt{[\Sigma_Y]_{jj}}$$

As you can see, the non-diagonal entries of $\Sigma_X, \Sigma_Y$ don't enter the calculations at all. So I'd think you cannot write this in terms of standard matrix multiplications.

However, if you have a $Diag$ function that keeps the diagonal entries but zeroes out all the non-diagonal entries, then you can write:

$$Cov(X,Y) = \sqrt{Diag(\Sigma_X)} \ P_{XY} \ \sqrt{Diag(\Sigma_Y)}$$

The multiplication is in that order because left-multiplying by a diagonal matrix scales the rows and right-multiplying scales the columns. (I know some typical software packages have a similar function although in the case of python np.diag you actually have to call it twice, so be careful!)