Question about Gaussians and Joint Distributions Question: Let $\mu \in \mathbb{R}^{m}$, and $\Sigma,\Sigma' \in \mathbb{R}^{m\times m}$. Let $X$ be an m-dimensional random vector with $X \sim \mathcal{N}
(\mu, \Sigma)$, and let $Y$ be a m-dimensional random vector such that $Y | X \sim \mathcal{N}(X, \Sigma')$. Derive the distribution and parameters for the joint distribution for the pair $(X,Y)$
Looking at the solution, it seems like first I have to establish the unconditional distribution of Y, and it's parameters. The solution suggests that I can construct $Y$ via $Y = X + Z$ where $Z \sim \mathcal{N}(0, \Sigma')$. My first question is, how come I am able to do this? What tells me that Y is a sum of two Gaussians? If I understand correctly, under closure properties of Gaussians, if $A$ and $B$ are both Gaussian, then $A+B$ is Gaussian, but I'm not sure how to deal with the conditional in this case. 
After this is established, the solutions provide that because $X, Y$ and $Y|X$ are all Gaussian, the joint distribution of $(X, Y)$ is also Gaussian. Why is this the case? It's not true that if $X$ and $Y$ are Gaussian, that $X*Y$ must also be, right? So how can I use the conditional to help me here?
Lastly, when the solution attempts to solve for the parameters, for the covariance, it will be in block form: $$\Sigma_{X,Y}=$$
\begin{bmatrix}
\Sigma_{X,X}  &\Sigma_{X,Y} \\
\Sigma_{Y,X}  &\Sigma_{Y,Y}
\end{bmatrix}
I'm not understanding how to calculate the two blocks between $X$ and $Y$.
I'm quite a beginner at this material, so more simple explanations appreciated. 
 A: First of all, what your solution should really say is that $Y \mid X$ is $$Y \mid X = X + (Z \mid X)$$
where $Z \mid X \sim \mathcal{N}(0, \boldsymbol\Sigma^{\prime})$. This wouldn't make sense otherwise. 
Why? Because
$$\mathbb{E}[Y \mid X] = \mathbb{E}[X+Z \mid X] = \mathbb{E}[X \mid X]+\mathbb{E}[Z \mid X] = X+\mathbf{0}=X$$
and $$\mathrm{Var}(Y \mid X) = \mathrm{Var}(X+(Z \mid X) \mid X)=\mathrm{Var}((Z \mid X)\mid X)=\boldsymbol\Sigma^{\prime}\text{.}$$
If this isn't satisfying enough for you, I would recommend going back to the definition of a multivariate normal random vector and then reconstructing $Y \mid X$ and $X$ from this, obtaining $(Z \mid X)$ as another multivariate normal random vector.
As for your second question, recall from probability that
$$f_{Y \mid X}(y \mid x) = \dfrac{f_{X, Y}(x, y)}{f_{X}(x)}\text{.}$$
Thus, you can derive the joint density 
$$f_{X, Y}(x, y)=f_{Y \mid X}(y \mid x)f_{X}(x)\text{.}$$
Plug in the appropriate densities above (i.e., $Y \mid X$ is multivariate normal, $X$ is multivariate normal), simplify it so that $f_{X, Y}$ looks like a multivariate normal, and find the appropriate covariance matrix using this.
