Joint trivariate normal density I now introduce the problem:
Let assume $\mathbf{z} = (z_1, z_2,z_3)$ be a trivariate normal variable. I want to find the covariance matrix of $\mathbf{z}$. 
I now that the density of $(z_1, z_2)$ is a bivariate normal   with mean vector $(\mu_1, \mu_2)$ and covariance matrix
$$
\left(
\begin{array}{cc}
\sigma_1^2 & \sigma_1\sigma_2 \rho \\
\sigma_1\sigma_2 \rho & \sigma_2^2 \\
\end{array} \right)
$$
and i now that $z_3 | z_1, z_2 = \beta_0+\beta_1z_1+\beta_2 z_2+\epsilon$ where $\epsilon \sim N(0, \sigma_3^2)$ and then the distribution of $z_3 | z_1, z_2$ is normal with mean $\beta_0+\beta_1z_1+\beta_2 z_2$ and variance  $\sigma_3^2$.
I want to find the distribution of $z_1, z_2, z_3$. I think that is a trivariate normal with mean $\mu_1, \mu_2, \beta_0$ and covariance matrix
$$
\left(
\begin{array}{ccc}
\sigma_1^2 & \sigma_1\sigma_2 \rho& \beta_1\sigma_1^2+\beta_2 \sigma_1\sigma_2 \rho  \\
\sigma_1\sigma_2 \rho & \sigma_2^2 &  \beta_1\sigma_1\sigma_2 \rho +\beta_2\sigma_2^2 \\
\beta_1\sigma_1^2+\beta_2 \sigma_1\sigma_2 \rho &  \beta_1\sigma_1\sigma_2 \rho +\beta_2\sigma_2^2  & \beta_1^2\sigma_1^2+\beta_2^2\sigma_2^2+\sigma_3^2
\end{array} \right)
$$
I think that this is right but with $\sigma_1^2=0.4$, $\sigma_2^2=1$, $\sigma_3^2=0.4$, $\rho=-0.1$, $\beta_1=-2$ and $\beta_2 = 3$ the matrix is not positive definite definite, i.e if i try to simulate a trivariate normal variable in R i get the following message:

Error in chol.default(V):
    the leading minor of order 3 is not positive definite

Where is the error?
 A: One can also use symmetry and linearity of the covariance operator, i.e.
$$\eqalign{\text{cov}(x+y,z)&=\text{cov}(x,z)+\text{cov}(y,z), \hspace{0.5em}\text{cov}(\alpha x,y) &= \alpha \,\text{cov}(x,y), \hspace{0.5em} \text{cov}(x,y) = \text{cov}(y,x)}$$
The missing element $\text{var}(z_3)$ can be found by
$$\eqalign{\text{var}(z_3)=\text{cov}(z_3,z_3)&=\text{cov}(\beta_0+\beta_1 z_1 + \beta_2 \,z_2 + \varepsilon,\beta_0+\beta_1 z_1 + \beta_2 \,z_2 + \varepsilon)\\
&=\beta_1^2 \text{cov}(z_1,z_1) + \beta_2^2  \text{cov}(z_2,z_2) + 2\beta_1 \beta_2 \text{cov}(z_1,z_2) + \text{cov}(\varepsilon,\varepsilon)\\
&= \beta_1^2 \sigma_1^2 + \beta_2^2\sigma_2^2 + 2\beta_1\beta_2\sigma_1\sigma_2\rho + \sigma_3^2}$$
where the covariances with $\beta_0$ and $\beta_1$ as arguments were dropped, because the covariance where one argument is constant yield zero. Furthermore $\text{cov}(\varepsilon,z_1)=0$ and $\text{cov}(\varepsilon,z_2)=0$
A: The Law of Total Covariance applied to $z_3$ asserts
$$\text{Var}(z_3) = \mathbb{E}(\text{Var}(z_3\ |\ (z_1,z_2)) + \text{Var}\left(\mathbb{E}(z_3\ |\ (z_1, z_2))\right)$$
whence, because $\mathbb{E}(\varepsilon)=0$ and $\text{Var}(\varepsilon)=\sigma_3^2$,
$$\eqalign{
\text{Var}(z_3) &= \mathbb{E}(\sigma_3^2) + \text{Var}(\mathbb{E}(\beta_0+\beta_1z_1+\beta_2z_2+\varepsilon\ |\ (z_1,z_2))) \\
&= \sigma_3^2 + \text{Var}(\beta_0 + \beta_1z_1 + \beta_2z_2) \\
&= \sigma_3^2 + \beta_1^2\sigma_1^2 + \beta_2^2\sigma_2^2 + 2\beta_1\beta_2\sigma_1\sigma_2\rho.
}$$
That is what belongs in the lower right entry of the covariance matrix.
