Consider i.i.d. random variables, $Z_1, \ldots ,Z_n \sim (0,1)$, written $\underline{Z}=\left( \begin{array}{ccc}
   Z_1 \\
  \vdots\\
  Z_n\\
 \end{array} \right)$ and let $\underline{Y}=A \underline{Z} + \underline{\mu}$ where $A$ is an invertible $n x n$ matrix and $\underline{\mu} \in \mathbb{R}^n$ is a vector, so $ Z= A^{-1}(Y-\underline{\mu})$.\\

Then the p.d.f. of $Y$ is given by $$f_{\underline{Y}}(\underline{y})= f_{\underline{Z}}(A^{-1}(\underline{y}- \underline{\mu})) \vert det(A^{-1}) \vert$$

But the joint p.d.f. of $\underline{Z}$ is the product of the p.d.f.'s of $Z_1, \ldots , Z_n$, so $f_{\underline{Z}}(\underline{z})= f(z_1) \cdot f(z_2) \cdot \ldots \cdot f(z_n)$ where

$$f(z_i) = \frac{1}{\sqrt{2 \pi}} e^{-\frac{z^2}{2}}$$

and hence

$$f_{\underline{Z}}(\underline{z}) = \prod_{i=1}^n \frac{1}{\sqrt{2 \pi}} e^{\frac{-z^2}{2}}$$

$$ = (\frac{1}{\sqrt{2 \pi}})^n e^{-\frac{1}{2} \sum_{i=1}^n z_i^2}$$

$$=\frac{1}{(2 \pi)^\frac{n}{2}} e^{-\frac{1}{2} \underline{z}'\underline{z}}$$

since

$$\sum_{i=1}^n z_i^2 = \Vert \underline{z} \Vert ^2 = \underline{z} \cdot \underline{z} = \underline{z}' \underline{z}$$

The joint p.d.f. of $\underline{Y}$ is therefore

$$f_{\underline{Y}}(\underline{y}) = f_{\underline{Z}}(A^{-1}(\underline{y} - \underline{\mu})) \vert det(A^{-1}) \vert $$

$$=\frac{1}{(2 \pi)^{\frac{n}{2}}} e^{-\frac{1}{2}(A^{-1}(\underline{y}-\underline{\mu}))'(A^{-1}(\underline{y}-\underline{\mu}))}\frac{1}{\vert det(A)\vert}$$

We can write $det(AA')=det(A)^2$ so $\vert det(A)\vert = \sqrt{det(AA')}$ and if we write $\Sigma=AA'$, then
$$\vert det(A) \vert = \vert \boldsymbol{\Sigma} \vert ^\frac{1}{2}$$


Also, note that
$$(A^{-1}(\underline{y}-\underline{\mu}))'(A^{-1}(\underline{y}-\underline{\mu})) = (\underline{y} - \underline{\mu})'(A^{-1})' A^{-1}(\underline{y} - \underline{\mu}) = (\underline{y} - \underline{\mu})' \boldsymbol{\Sigma}^{-1}(\underline{y}-\underline{\mu})$$

We can now write

$$f_{\underline{Y}}(\underline{y}) = \frac{1}{(2 \pi)^\frac{n}{2} \vert \boldsymbol{\Sigma} \vert ^{\frac{1}{2}}} e^{-\frac{1}{2} (\underline{y}-\underline{\mu}) \boldsymbol{\Sigma}^{-1} (\underline{y}-\underline{\mu})}$$

This is the density of the multivariate normal distribution.

Note that
$$E[\underline{Y}] = \mu$$
$$V[\underline{Y}] = V[A\underline{Z}] = AV[\underline{Z}]A' = AIA' = \boldsymbol{\Sigma}$$

Notation: $\underline{Y}\sim n(\underline{\mu}, \boldsymbol{\Sigma})$