I try to understand the joint distribution of two random variables following a normal distribution. Wikipedia gives the result but I don't catch how do they get the equation:
The probability density function of a vector $[X\;Y]$ is :
$ f(x,y) = \frac{1}{2 \pi \sigma_X \sigma_Y \sqrt{1-\rho^2}} \exp\left( -\frac{1}{2(1-\rho^2)}\left[ \frac{(x-\mu_X)^2}{\sigma_X^2} + \frac{(y-\mu_Y)^2}{\sigma_Y^2} - \frac{2\rho(x-\mu_X)(y-\mu_Y)}{\sigma_X \sigma_Y} \right] \right)$
Let's $X$ and $Y$ be two random variable following a normal distribution with respectively $\mu_x$ and $\mu_y$ as mean and $\sigma_x$ and $\sigma_y$ as standard deviation:
$P(X=x) = \frac{1}{\sqrt{2\pi}\sigma_x}exp(-\frac{1}{2}(\frac{x-\mu_x}{\sigma_x})^2)$
and:
$P(Y=y) = \frac{1}{\sqrt{2\pi}\sigma_y}exp(-\frac{1}{2}(\frac{y-\mu_y}{\sigma_y})^2)$
If we suppose that $X$ and $Y$ are uncorrelated I understand that we can derive the joint probability distribution of $X$ and $Y$ ($P(X=x \; and \; Y=y)$) by multiplying the two distributions:
$P(X=x \; and \; Y=y) = \frac{1}{\sqrt{2\pi}\sigma_x}exp(-\frac{1}{2}(\frac{x-\mu_x}{\sigma_x})^2) \times \frac{1}{\sqrt{2\pi}\sigma_y}exp(-\frac{1}{2}(\frac{y-\mu_y}{\sigma_y})^2)$
$= \frac{1}{2\pi\sigma_x\sigma_y}exp(-\frac{1}{2}(\frac{x-\mu_x}{\sigma_x})^2-\frac{1}{2}(\frac{y-\mu_y}{\sigma_y})^2)$
$= \frac{1}{2\pi\sigma_x\sigma_y}exp(-\frac{1}{2}((\frac{x-\mu_x}{\sigma_x})^2+(\frac{y-\mu_y}{\sigma_y})^2)))$
But how are we suppose to manage the fact that $X$ and $Y$ could be correlated ?
Well, you can't manage the fact that $X$ and $Y$ could be correlated. Or, rather, simply knowing the marginals of $X$ and $Y$ and their correlation is not enough to know the joint distribution $(X,Y)$, even if we know that both $X$ and $Y$ are normal.
To see this, let $S\sim Ber(p)$ for $p\in [0,1]$ and assume that $X\sim \mathcal{N}(0,1)$ and define $Y=(-1)^SX$. Assuming $S$ and $X$ are independent, then $\mathbb{P}(Y\leq t)=\mathbb{P}(S=0)\mathbb{P}(X\leq t)+\mathbb{P}(S=1)\mathbb{P}(X\geq -t)=\mathbb{P}(X\leq t),$ since $\mathcal{N}(0,1)$ is symmetric. Hence, $Y\sim \mathcal{N}(0,1)$. Furthermore, $Cov(X,Y)=1-2p$ by direct computation, so doing this experiment, we can get $Cov(X,Y)$ to be anything we want - for instance, they could be uncorrelated without being independent.
However, the distribution of $(X,Y)$ does not have density (with respect to the two-dimensional Lesbegue measure), since for $A=\{(t,t)|t\in \mathbb{R}\}\cup \{(-t,t)|t\in \mathbb{R}\},$ we have $\mathbb{P}((X,Y)\in A)=1$, but $A$ is a set of measure $0$.
So what's wrong? Well, we can't just pretend we know that $X$ and $Y$ are normal and have some covariance. It is absolutely vital that we know that $(X,Y)$ is a (regular) normal vector. And this is just defined as a variable such that there exists iid normal variables $Z_1,Z_2$ and an invertible matrix $B$ and $\mu\in \mathbb{R}^2$ such that $(X,Y)=B(Z_1,Z_2)+\mu$.
In this case, we can use the multidimensional change of variables formula (https://en.wikipedia.org/wiki/Integration_by_substitution#Substitution_for_multiple_variables) to deduce that $(X,Y)$ has density $$ \frac{1}{2\pi \det(BB^T)}\exp\left(-\frac{1}{2}(x-\mu)^T (B B^T)^{-1}(x-y)\right), $$ which agrees with what you have above.