Covariance matrix for a linear combination of correlated Gaussian random variables Supposing $X$ and $Y$ are random variables with a joint bivariate normal distribution and covariance matrix $\Sigma_{XY}$. Consider the following linear combination for constants $A$, $B$ and $C$:
$$Z = AX + BY + C$$
This Cross Validated post states that $Y$ will be normally distributed according to:
$$
AX+BY+C \sim \mathcal{N}\left[ 
\left(\begin{matrix}A& B \end{matrix}\right)
\left(\begin{matrix}\mu_X\\\mu_Y\end{matrix}\right) + C, 
\left(\begin{matrix}A & B \end{matrix}\right)\Sigma_{X,Y} \left(\begin{matrix}A^T \\ B^T \end{matrix}\right)\right]
$$
However I would like to know how to compute the covariance matrix $\Sigma_{X,Y,Z}$. How can I calculate the correlation between the new variable $Z$ and the old variables $X$ and $Y$ i.e. $corr(Z,X)$ and $corr(Z,Y)$?
 A: If $X$ and $Y$ are correlated (univariate) normal random variables and $Z = AX+BY+C$,  then the linearity of expectation and the bilinearity of the covariance function gives us that
\begin{align}
E[Z] &= AE[X] + BE[Y] + C,\tag{1}\\
\operatorname{cov}(Z,X) &= \operatorname{cov}(AX+BY+C,X) = A\operatorname{var}(X) + B\operatorname{cov}(Y,X)\\
\operatorname{cov}(Z,Y) &= \operatorname{cov}(AX+BY+C,Y) = B\operatorname{var}(Y) + A\operatorname{cov}(X,Y)\\
\operatorname{var}(Z) &= \operatorname{var}(AX+BY+C) \quad = A^2\operatorname{var}(X) + B^2\operatorname{var}(Y) + 2AB \operatorname{cov}(X,Y), \tag{2}\\
\end{align}
but it is not necessarily true that $Z$ is a normal (a.k.a Gaussian)
random variable.  That $X$ and $Y$ are jointly normal random variables
is sufficient to assert that $Z = AX+BY+C$ is a normal random variable.
Note that $X$ and $Y$ are not required to be independent; they can be
correlated as long as they are jointly normal. For examples of
normal random variables $X$ and $Y$ that are not jointly normal
and yet their sum $X+Y$ is normal, see the answers to 
Is joint normality a
necessary condition for the sum of normal random variables to 
be normal?. As pointed out at the end of my own answer there, joint normality means that
all linear combinations $aX+bY$ are normal, whereas in the special
case being discussed there, only one linear combination
$X+Y$ of non-jointly normal random variables is proven to be normal; most other linear combinations are not normal.
More generally, if $X$ and $Y$ are (column) $n$-vector random variables with $n\times n$ covariance matrices $\Sigma_{X,X}$, $\Sigma_{Y,Y}$,
and $n\times n$ crosscovariance matrix $\Sigma_{X,Y}$,
$A$ and $B$ are $m\times n$ nonrandom matrices, and $Z$ and $C$ (column)
$m$-vectors, then it is indeed true that
\begin{align} E[Z] &= AE[X] + BE[Y] + C &\quad \scriptstyle{\text{compare with } (1)}\\
\Sigma_{Z,Z} &= A\Sigma_{X,X}A^T + B\Sigma_{Y,Y}B^T +2A\Sigma_{X,Y}B^T
&\quad \scriptstyle{\text{compare with } (2)}\\
\end{align}
but, as in the univariate case, it is not necessarily true that
$Z$ is a normal vector (in the sense that the $m$ components $Z_i$
are jointly normal random variables). Once again, joint
normality of $(X_1, X_2, \ldots, X_n, Y_1, Y_2, \ldots, Y_n)$
suffices to allow the assertion that $Z$ is a normal random vector.
A: Covariance matrix $\Sigma_{xy}$ can be written:
$$ \Sigma_{xy} = \left[ \begin{array}{cc} c_{xx} & c_{xy} \\ c_{yx} & c_{yy}  \end{array} \right] $$
Where $c_{ab}$ denotes the covariance between $a$ and $b$. (Note symmetric so $c_{ab} = c_{ba}$.)
The covariance matrix for $X$,$Y$, and $Z$ would be:
$$ \Sigma_{xyz} = \left[ \begin{array}{ccc} c_{xx} & c_{xy} & c_{xz} \\ c_{yx} & c_{yy} & c_{yz} \\ c_{zx} & c_{zy} & c_{zz} \\  \end{array} \right] $$
You need to find the additional terms:
$$\begin{align*}
c_{zx} &= E[(Z - E[z])(X - E[x])] = \quad ?\\
c_{zy} &= E[(Z - E[z])(Y - E[y])] = \quad ?\\
c_{zz} &= E[(Z - E[z])^2]  = \quad ?\\
\end{align*}$$
