Multivariate Normal Difference Distribution

Solution 1:

If two multivariate normal random variables $X\sim\mathcal N(\mu_X,\Sigma_X)$ and $Y\sim\mathcal N(\mu_Y,\Sigma_Y)$ are independent (and of the same dimension), then their sum is still normal, and you can sum the mean and variance directly: $$X+Y\sim\mathcal N(\mu_X+\mu_Y,\Sigma_X+\Sigma_Y).$$

The easiest way to see this is with the expression of the characteristic function, as can be found on the Wikipedia article for example: $$\phi_X(u)=\exp(iu'\mu_X-u'\Sigma_Xu).$$

Since the characteristic function of the sum of two independent random variables is the product of their respective characteristic functions, then \begin{eqnarray} \phi_{X+Y}(u) &=& \exp(iu'\mu_X-u'\Sigma_Xu)\exp(iu'\mu_Y-u'\Sigma_Yu)\\ &=& \exp(iu'\mu_X-u'\Sigma_Xu+iu'\mu_Y-u'\Sigma_Yu)\\ &=& \exp(iu'(\mu_X+\mu_Y)-u'(\Sigma_X+\Sigma_Y)u). \end{eqnarray}

Now back to your original question, it suffices to see that $-Y$ has the same variance as $Y$, but with opposite mean to conclude: $$-Y\sim \mathcal N(-\mu_Y,\Sigma_Y),$$ so indeed $$X-Y\sim\mathcal N(\mu_X-\mu_Y,\Sigma_X+\Sigma_Y).$$