From what I understand drifted Brownian motion as an SDE in one dimension has a normal distribution. In two dimensions the SDE
\begin{align} dX_t = \mu_x X_t + \sigma_x dW^x_t \\ dY_t = \mu_y Y_t + \sigma_y dW^y_t, \end{align} would follow the bivariate normal distribution, from the definition of independence $p(x,y) = p(x)p(y)$.
My question is how would we derive the density function that would arise if we had correlated noise/Wiener processes $W^x$ and $W^y$?
Can it still be the bivariate normal distribution when the processes are not independent? If not, how would you derive/explain this?
For an $n$-dimensional Ornstein-Uhlenbeck process \begin{equation} dX_{t}=AX_{t}dt+\sqrt{C}dB_{t}, \quad X_{0}=x \end{equation} with $A\in\mathbb{R}^{n\times n}$ (which is $\operatorname{diag}(\mu_{x},\mu_{y})$ in your case) and correlation matrix $\sqrt{C}\in\mathbb{R}^{n\times n}$ solutions are given by \begin{equation} X_{t}=e^{tA}x+\int_{0}^{t}e^{(t-s)A}\sqrt{C}dB_{s} \end{equation} from which it follows by some linear algebra and standard properties of Ito-Integrals that $X_{t}$ is $\mathcal{N}(e^{tA}x,Q_{t})$ normally distributed, where \begin{equation} Q_{t}=\int_{0}^{t}e^{sA}Ce^{sA^{*}}ds. \end{equation}