Suppose $X$ is a $p$-dimensional multivariate normal random variable with conditional distribution
$$f_X (x\mid Z=z) = \mathcal{N}(\mu + Az, I\sigma^2),$$
Further suppose that
$$Z\sim \mathcal{N}(0, I),$$
where $Z$ is a $q$-dimensional multivariate normal random variable. (So in fact, $A$ is a $p\times q$ matrix.)
What is the joint distribution of $X$ and $Z$? In other words, what is the pdf
$$f_{X_1,\cdots,X_p,Z_1,\cdots,Z_q}(x_1,\cdots,x_p,z_1,\cdots,z_q)?$$
We can assume all the independence that we want (i.e. all $x_i$'s and $z_j$'s are independent)
I spent quite some time googling but all the examples I came across seem to implicitly assume that $X$ and $Y$ have the same dimension. For example, in this Matrix Cookbook page 41, it talks about "product of Gaussian density" but assumes that two input vectors need to have same dimensions.
Best Answer
$$ X \mid Z=z \sim \mathcal{N}(\mu + Az, \sigma^2 I) $$ $$ (X- AZ) \mid (Z=z) \sim\mathcal{N}(\mu,\sigma^2 I) $$ Since the conditional distribution of $X-AZ$ given $Z$ does not depend on $Z,$ we must conclude that (1) $X-AZ$ and $Z$ are independent, and (2) the marginal (or "unconditional") distribution of $X-AZ$ is that same distribution.
Thus we have $$ X-AZ\sim\mathcal N(\mu,\sigma^2 I) $$ and $$ (X-Az) \text{ and } Z \text{ are independent.} $$ Thus $$ \begin{bmatrix} X-AZ \\ Z \end{bmatrix} \sim \mathcal N \left( \begin{bmatrix} \mu \\ 0 \end{bmatrix}, \begin{bmatrix} \sigma^2 I & 0 \\ 0 & I \end{bmatrix} \right). $$ $$ \begin{bmatrix} X \\ Z \end{bmatrix} = \begin{bmatrix} I & A \\ 0 & I \end{bmatrix} \begin{bmatrix} X-AZ \\ Z \end{bmatrix} $$
Therefore $\begin{bmatrix} X \\ Z \end{bmatrix}$ is multivariate normal and $$ \operatorname{E}\begin{bmatrix} X \\ Z \end{bmatrix} = \begin{bmatrix} I & A \\ 0 & I \end{bmatrix} \operatorname{E} \begin{bmatrix} X-AZ \\ Z \end{bmatrix} = \begin{bmatrix} I & A \\ 0 & I \end{bmatrix} \begin{bmatrix} \mu \\ 0 \end{bmatrix} = \begin{bmatrix} \mu \\ 0 \end{bmatrix}, $$ and $$ \operatorname{var} \begin{bmatrix} X \\ Z \end{bmatrix} = \begin{bmatrix} I & A \\ 0 & I \end{bmatrix} \begin{bmatrix} \sigma^2 I & 0 \\ 0 & I \end{bmatrix} \begin{bmatrix} I & 0 \\ A^T & I \end{bmatrix} = \begin{bmatrix} \sigma^2 I + A A^T & A \\ A^T & I \end{bmatrix}. $$