Prove sum of unindependent Gaussian random variable is still Gaussian I know how to prove the statement when two r.v. are independent. But how to prove it when the r.v. are correlated?
Edit:or counter example?
 A: In some cases in which they are correlated it is not true that their sum has a normal (or "Gaussian") distribution.  For example, suppose $X\sim N(0,1)$ and let $Y = \pm X$, where $+$ or $-$ is chosen randomly, and you pick "$+$" with probability $2/3$ and "$-$" with probability $1/3$, independently of the value of $X$.  Then it's not hard to show that $Y\sim N(0,1)$, and $X$ is positively correlated with $Y$.  But $\Pr(X+Y=0) = 1/3$, so $X+Y$ cannot be normally distributed.
If $X$ and $Y$ are jointly normally distributed, then it is indeed true that $X+Y$ is normally distributed regardless of the correlation between $X$ and $Y$.  That's a trivial conseequence of a statement often taken to be the definition of joint normality: $(X,Y)$ has a bivariate normal distribution (so $X$ and $Y$ are "jointly normally distributed") precisely if the distribution of the pair $(X,Y)$ is such that every linear combination $aX+bY$, where $a,b$ are not random, is normally distributed.
Another definition of joint normality says there are standard normal random variables $Z_1,\ldots,Z_k$ such that $X$ is some linear combination of them, and $Y$ is some other linear combination of them.  Then $X+Y$ is normally distributed because it's a linear combination of indendent normals $Z_1,\ldots,Z_k$.
Another characterization says the joint density is
$$
\text{constant}\cdot\exp\left( \frac{-1}{2} [x-\mu,y-\nu]\begin{bmatrix} a & b \\ b & c \end{bmatrix}^{-1}\begin{bmatrix} x-\mu \\ y-\nu \end{bmatrix} \right)
$$
where the matrix is a positive-definite.  This is not quite equivalent to the foregoing to characterizations because it does not accomodate the case where the correlation is either $1$ or $-1$, or, in the case of a random vector with more than two components, the case where the matrix that is the variance is singular.  Proving $X+Y$ is normally distributed in this case is something I would do by showing that if $\mathbf X\in\mathbb R^n$ and $\mathbf X\sim N_n(\mu,A)$, and $M\in\mathbb R^{k\times n}$ is a constant (i.e. not random) matrix, then $M\mathbf X \sim N_k(M\mu,MAM^\top)$.  Maybe I'll add more about that later.
A: In general, for jointly continuous random variables $X$ and $Y$ with
joint pdf $f_{X,Y}(u,v)$,$$f_{X+Y}(\alpha) = \int_{-\infty}^\infty f_{X,Y}(u, \alpha-u)\,\mathrm du$$
which reduces to the convolution of the marginal densities $f_X$ and $f_Y$
when $X$ and $Y$ are independent random variables.  For the special case
when $X$ and $Y$ are jointly Gaussian random variables, the integral
above can be computed by "completing the square" in the exponent so that you
can pull out some things not dependent on $u$ outside the integral and leave
a Gaussian density inside the integral, which then integrates to $1$.  I assume
that you know what a bivariate joint Gaussian density is.
If $X \sim N(\mu_X, \sigma_X^2)$ and $Y \sim N(\mu_Y, \sigma_Y^2)$
are jointly Gaussian random variables with correlation coefficient $\rho$,
then $$X+Y \sim N(\mu_X+\mu_Y,\sigma_X^2+\sigma_Y^2+2\rho\sigma_X\sigma_Y)$$
