Constructing a bivariate normal from three univariate normals I'm trying to construct correlated bivariate normal random variables from three univariate normal random variables. I realize there is a formula for constructing a bivariate normal random variable from two univariate random normal variables, but I have reasons for wanting to adjust two previously sampled variables by a third in order to give them a correlation $\rho$. 
Based on the approaches for constructing bivariate normals from two univariate normals, I came up with the following approach and need help verifying its correctness. 
First, imagine that we have three univariate normal variables. For simplicity here, we just assume they all have $\sigma=1$.
$$ X_0 \sim Normal(0, 1)$$
$$ Y_0 \sim Normal(0, 1)$$
$$ Z \sim Normal(0, 1)$$
Given these three univariate random variables, I construct two new random variables using the following linear combinations:
$$ X = |\rho| * Z + \sqrt{1-\rho^2} * X_0$$
$$ Y = \rho * Z + \sqrt{1-\rho^2} * Y_0$$
where $\rho \in [-1, 1]$ represents the correlation coefficient between the two univariate normals.
Can someone help me formally verify that $X$ and $Y$ are now correlated random variables with correlation $\rho$? 
I've convinced my self through empirical simulation. Here are plots of values sampled from $X$ and $Y$ for the cases where $\rho=0$, $\rho=1$, and $\rho=-1$. 
Plot of X vs. Y for rho of 0
Plot of X vs. Y for rho of 1
Plot of X vs. Y for rho of -1
These plots are exactly as I would expect, but it would be nice to have a formal proof based on my construction. Thanks in advance!
 A: I presume $X_0, Y_0$ and $Z$ are assumed to be independent? Well, if we define
$$ \mathbf U := \begin{bmatrix} X_0 \\ Y_0 \\Z \end{bmatrix}, \ \ \ \ \ \ \mathbf V := \begin{bmatrix} X \\ Y\end{bmatrix},$$
then
$$ \mathbf V = \mathbf M \mathbf U,$$
where
$$ \mathbf M := \begin{bmatrix} \sqrt{1- \rho^2} & 0 & | \rho| \\ 0 & \sqrt{1- \rho^2} & \rho \end{bmatrix}$$
Since $X_0, Y_0$ and $Z$ are independent with unit variance, the covariance matrix for $\mathbf U$ is the $3 \times 3$ identity matrix:
$$\mathbb E[\mathbf U \mathbf U^T] = \mathbf I.$$
Hence the covariance matrix for $\mathbf V$ is given by
$$ \mathbb E [\mathbf V \mathbf V^T ] = \mathbf M \ \mathbb E[\mathbf U \mathbf U^T ] \ \mathbf M^T = \mathbf M\mathbf M^T=\begin{bmatrix} 1 & \rho | \rho| \\ \rho | \rho| & 1\end{bmatrix}$$
Thus the correlation coefficient between $X$ and $Y$ is
$$ \rho_{X,Y} = \frac{\mathbb E[XY]}{\sqrt{\mathbb E[X^2] \ \mathbb E[{Y^2]}}} = \frac{\rho| \rho|}{\sqrt{1 \times 1}} = \rho|\rho|.$$
[I used the fact that $\mathbb E[X] = \mathbb E[Y] = 0$ here.]
So I'm afraid the correlation coefficient is not $\rho$. It's ${\rm sign}(\rho) \times  |\rho|^2$.
But that's easily fixed. If you redefine $\mathbf M$ as $$ \mathbf M := \begin{bmatrix} \sqrt{1- |\rho |} & 0 &  {\rm sign}(\rho)\sqrt{ |\rho |} \\ 0 & \sqrt{1- |\rho | } & \sqrt{|\rho |} \end{bmatrix},$$
then it should work out.
Finally, I'll point out that the $X$ and $Y$ that you've constructed are guaranteed to be Gaussian, since $\mathbf V$ is related to the Gaussian vector $\mathbf U$ by a linear transformation. They also both have unit variance, as you can see from the diagonal elements of $\mathbb E[\mathbf V \mathbf V^T]$.
