What is the analytic expression for PDF of joint distribution of two Gaussian random vectors?

Question

I know that if $X$ and $Y$ are random variables with respective PDFs,

$$ f_X(x) = \frac{1}{\sqrt{2\pi\sigma_x^2}}\exp\left\{-\frac{\left(x-\mu_x\right)^2}{2\sigma_x^2}\right\} \\ f_Y(y) = \frac{1}{\sqrt{2\pi\sigma_y^2}}\exp\left\{-\frac{\left(y-\mu_y\right)^2}{2\sigma_y^2}\right\} $$

Then their joint PDF is written as

$$ f_{XY}(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{2\rho(x-\mu_x)(y-\mu_y)}{\sigma_x \sigma_y} + \frac{(y-\mu_y)^2}{\sigma_y^2} \right] \right) $$

But when $\mathbf{x}$ and $\mathbf{y}$ are random vectors with PDFs

$$ f_{\mathbf x}(x_1,\ldots,x_k)\, = \frac{1}{(2\pi)^{k/2}|\boldsymbol\Sigma_x|^{1/2}} \exp\left(-\frac{1}{2}({\mathbf x}-{\boldsymbol\mu_x})^T{\boldsymbol\Sigma_x}^{-1}({\mathbf x}-{\boldsymbol\mu_x}) \right) \\ f_{\mathbf y}(y_1,\ldots,y_k)\, = \frac{1}{(2\pi)^{k/2}|\boldsymbol\Sigma_y|^{1/2}} \exp\left(-\frac{1}{2}({\mathbf y}-{\boldsymbol\mu_y})^T{\boldsymbol\Sigma_y}^{-1}({\mathbf y}-{\boldsymbol\mu_y}) \right) $$

How do you express their joint PDF?

$$ f_{\mathbf xy}(\mathbf {x,y})\, = \, ? $$

Answer 1

In the case where you only assume that $\mathbf{X}$ and $\mathbf{Y}$ are marginally Gaussian, you can't say much about the joint density of $(\mathbf{X},\mathbf{Y})$, and you certainly can't conclude that the joint density is a Gaussian density. In the answer below I've added the additional assumption that the joint distribution is indeed Gaussian.

Assume that $\mathbf{X}=(X_1,\ldots,X_k)\sim\mathcal{N}_k(\boldsymbol{\mu}_x,\mathbf{\Sigma}_x)$ and $\mathbf{Y}=(Y_1,\ldots,Y_k)\sim\mathcal{N}_k(\boldsymbol{\mu}_y,\mathbf{\Sigma}_y)$ and that $\mathbf{Z}=(\mathbf{X},\mathbf{Y})$ follows a Gaussian distribution, then the joint density is given by $$ f_\mathbf{Z}(\mathbf{z})= \frac{1}{(2\pi)^{2k/2}|\boldsymbol\Sigma_z|^{1/2}} \exp\left(-\frac{1}{2}({\mathbf z}-{\boldsymbol\mu_z})^T{\boldsymbol\Sigma_z}^{-1}({\mathbf z}-{\boldsymbol\mu_z}) \right),\quad\mathbf{z}\in\mathbb{R}^{2k}, $$ where $$ \boldsymbol{\mu_z}=(\boldsymbol{\mu}_x,\boldsymbol{\mu}_y) $$ and $$ \mathbf{\Sigma}_z=\begin{bmatrix} \mathbf{\Sigma}_x & \mathbf{\Sigma}_{xy} \\ \mathbf{\Sigma}_{yx} & \mathbf{\Sigma}_y \end{bmatrix} $$ written in block form. Here $\mathbf{\Sigma}_{xy}$ is the $k\times k$ matrix whose $(i,j)$th entry is $$ \mathbf{\Sigma}_{xy}^{i,j}=\mathrm{Cov}(X_i,Y_j). $$

Note that this also holds when $\mathbf{X}$ is $k$-dimensional and $\mathbf{Y}$ is $n$-dimensional for $n\neq k$.