What is the distribution of the sum of non i.i.d. gaussian variates? If $X$ is distributed $N(\mu_X, \sigma^2_X)$,
$Y$ is distributed $N(\mu_Y, \sigma^2_Y)$
and $Z = X + Y$, I know that $Z$ is distributed $N(\mu_X + \mu_Y, \sigma^2_X + \sigma^2_Y)$ if X and Y are independent.
But what would happen if X and Y were not independent, i.e.
$(X, Y) \approx N\big(
(\begin{smallmatrix}
\mu_X\\\mu_Y
\end{smallmatrix})
,
(\begin{smallmatrix}
\sigma^2_X && \sigma_{X,Y}\\
\sigma_{X,Y} && \sigma^2_Y
\end{smallmatrix})
\big)
$
Would this affect how the sum $Z$ is distributed?
 A: @dilip's answer is sufficient, but I just thought I'd add some details on how you get to the result.  We can use the method of characteristic functions.  For any $d$-dimensional multivariate normal distribution $X\sim N_{d}(\mu,\Sigma)$ where $\mu=(\mu_1,\dots,\mu_d)^T$ and $\Sigma_{jk}=cov(X_j,X_k)\;\;j,k=1,\dots,d$, the characteristic function is given by:
$$\varphi_{X}({\bf{t}})=E\left[\exp(i{\bf{t}}^TX)\right]=\exp\left(i{\bf{t}}^T\mu-\frac{1}{2}{\bf{t}}^T\Sigma{\bf{t}}\right)$$
$$=\exp\left(i\sum_{j=1}^{d}t_j\mu_j-\frac{1}{2}\sum_{j=1}^{d}\sum_{k=1}^{d}t_jt_k\Sigma_{jk}\right)$$
For a one-dimensional normal variable $Y\sim N_1(\mu_Y,\sigma_Y^2)$ we get:
$$\varphi_Y(t)=\exp\left(it\mu_Y-\frac{1}{2}t^2\sigma_Y^2\right)$$
Now, suppose we define a new random variable $Z={\bf{a}}^TX=\sum_{j=1}^{d}a_jX_j$.  For your case, we have $d=2$ and $a_1=a_2=1$.  The characteristic function for $Z$ is the basically the same as that for $X$.
$$\varphi_{Z}(t)=E\left[\exp(itZ)\right]=E\left[\exp(it{\bf{a}}^TX)\right]=\varphi_{X}(t{\bf{a}})$$
$$=\exp\left(it\sum_{j=1}^{d}a_j\mu_j-\frac{1}{2}t^2\sum_{j=1}^{d}\sum_{k=1}^{d}a_ja_k\Sigma_{jk}\right)$$
If we compare this characteristic function with the characteristic function $\varphi_Y(t)$ we see that they are the same, but with $\mu_Y$ being replaced by $\mu_Z=\sum_{j=1}^{d}a_j\mu_j$ and with $\sigma_Y^2$ being replaced by $\sigma^2_Z=\sum_{j=1}^{d}\sum_{k=1}^{d}a_ja_k\Sigma_{jk}$.  Hence because the characteristic function of $Z$ is equivalent to the characteristic function of $Y$, the distributions must also be equal.  Hence $Z$ is normally distributed.  We can simplify the expression for the variance by noting that $\Sigma_{jk}=\Sigma_{kj}$ and we get:
$$\sigma^2_Z=\sum_{j=1}^{d}a_j^2\Sigma_{jj}+2\sum_{j=2}^{d}\sum_{k=1}^{j-1}a_ja_k\Sigma_{jk}$$
This is also the general formula for the variance of a linear combination of any set of random variables, independent or not, normal or not, where $\Sigma_{jj}=var(X_j)$ and $\Sigma_{jk}=cov(X_j,X_k)$.  Now if we specialise to $d=2$ and $a_1=a_2=1$, the above formula becomes:
$$\sigma^2_Z=\sum_{j=1}^{2}(1)^2\Sigma_{jj}+2\sum_{j=2}^{2}\sum_{k=1}^{j-1}(1)(1)\Sigma_{jk}=\Sigma_{11}+\Sigma_{22}+2\Sigma_{21}$$
A: See my comment on probabilityislogic's answer to this question.  Here,
$$
\begin{align*}
X + Y &\sim N(\mu_X + \mu_Y,\; \sigma_X^2 + \sigma_Y^2 + 2\sigma_{X,Y})\\ 
aX + bY &\sim N(a\mu_X + b\mu_Y,\; a^2\sigma_X^2 + b^2\sigma_Y^2 + 2ab\sigma_{X,Y})
\end{align*}
$$ 
where $\sigma_{X,Y}$ is the covariance of $X$ and $Y$.  Nobody writes the off-diagonal entries in the covariance matrix as $\sigma_{xy}^2$ as you have 
done. The off-diagonal entries are covariances which 
can be negative.
