Posterior Predictive Distribution - Normal model with unknown mean, and known std dev I do not understand how one can prove that in a normal model, where we know the variance, but not the mean, the posterior predictive distribution is also normal...
Below is cut of a page from Gelman's Data Bayesian Analysis. I don't understand how the authors can state that the posterior predictive dist. will be a Normal one, i.e., why can they conclude that it's a joint normal posterior distribution and then conclude that the marginal is also normal?

Any help would be appreciated.
 A: if $X \sim N(\mu_X, \sigma_X^2)$ and $Y \sim N(\mu_Y, \sigma_Y^2)$ and $X$ and $Y$ are independent then 
$$
f_Z(z) = \int_{-\infty}^\infty f_Y(z-x) f_X(x) dx
$$
is a normal density with mean $\mu_X+\mu_Y$ and covariance $\sigma^2_X+\sigma_Y^2$ (that is actually the density of $Z=X+Y$). Simple proof at https://en.wikipedia.org/wiki/Sum_of_normally_distributed_random_variables#Proof_using_convolutions which states your result.
A: If the likelihood is normal, i.e. $y_i|\theta \sim N(\theta,\sigma_x^2)$ and the prior is normal, i.e. $\theta_{prior} \sim N(\theta_0, \sigma_{\theta}^2) $, then the posterior is also normal, $\theta \sim N(\mu_1, \tau_1^2)$. You can follow a proof here.
Now, the predictive distribution will also be normal:
$P(y'|y) = \int_{\theta}P(y', \theta |y)d\theta = \int_{\theta}P(y'|y,\theta)P(\theta|y)d\theta = \int_{\theta}P(y'|\theta)P(\theta|y)d\theta$
The first step is just treating the probability as a marginal probability over the parameter space. The 2nd step is just the conditional probability. The 3rd step is the assumption that once we know $\theta$, it's the only thing that influences the new observation, and any old observation is useless. What you end up with is the integral over all $\theta$ of the likelihood times the posterior.
You could calculate these integrals and see that you end up with a normal distribution. Or you can use a little trick, using convolutions (that peuhp hints of):
$\int_{\theta}P(y' = y_0|\theta)P(\theta|y)d\theta = \int_{\theta}P(y' = y_0 - \theta+\theta|\theta)P(\theta|y)d\theta = \int_{\theta}P(y' = y_1 + \theta|\theta)P(\theta|y)d\theta = \int_{\theta}P(y'-\theta = y_1|\theta)P(\theta|y)d\theta = \int_{\theta}P(y'-\theta|\theta)P(\theta|y)d\theta$
Here $y_0$ is the value that the random variable is equal to. I used variable change and called $y_1 = y_0-\theta$.
Note that what we end up with is a convolution of two functions which belong to two distributions: $\theta$ which we already know; and $y-\theta |\theta \sim N(0, \sigma_x^2)$. I.e. what we get is the sum of these two normal distributions, which is equal to:
$$ y'|y \sim N(\mu_1, \sigma_x^2 + \tau_1^2)
$$
A: I'm not sure the answers given so far reflect the author's reasoning, so I'll give my take on this question.
First off, I don't think the authors are using the approach of recognizing the integral as proportional to the pdf of the sum of two normal random variables. If they did that, there would be no need for them to derive the mean and variance using the law of iterated expectation and the law of total variance, as they do a bit further down the page.  But if one wants to use that approach, it does not seem like we can use the case for the sum of two independent normal random variables, which is what the two current answers seem to propose. As far as I can see, one would define a random variable $U$ such that $\tilde Y = \Theta + U$, and then use the fact that (see for example page 3 in this document): $$f_{\tilde Y}(\tilde y) = \int_{-\infty}^\infty f_{U|\Theta}(\tilde{y}-\theta)f_\Theta(\theta)d\theta \quad (*)$$
(the conditioning on $y$ done implicitly). The integral $(*)$ is what one gets when calculating the distribution of the sum of the two random variables $U$ and $\Theta$ where $U = \tilde{Y} - \Theta$, which given $\Theta$ has a distribution $\sim N(0, \sigma_{\tilde{Y}}^2)$ $\Big[$since $\tilde{Y}|\Theta \sim N(\Theta, \sigma_{\tilde{Y}}^2)$ $\Big]$. Note that $U$ and $\Theta$ are not independent.
Showing that this implies $\tilde Y \sim N(\mu_{\Theta}, \sigma_{\tilde{Y}}^2 + \sigma_{\Theta}^2)$ is done in the same way as in the case when $U$ and $\Theta$ are independent, already mentioned by the other answers, since the resulting integral $(*)$ still has an integrand which is the product of the same type of normal densities.

But the question asked was

...why can they conclude that it's a joint normal posterior distribution and then conclude that the marginal is also normal?

The second part of the question, how the authors can conclude that the marginal is also normal, is a consequence of a standard fact about the multivariate normal distribution: if $X$ is a random vector that follows the multivariate normal distribution, i.e. $X\sim N(\mu, \Sigma)$, then any subset vector $X_s$ follows the multivariate normal distribution, $X_s \sim N(\mu_s, \Sigma_s)$, where $\mu_s$ and $\Sigma_s$ drops the rows/columns not part of the subset. This is stated here and also on the wiki page.
For the first part of the question, how the authors can conclude that $\tilde{y}$ and $\theta$ have a joint normal distribution, the justification given by the authors from the quoted passage is that

The product in the integrand is the exponential of a quadratic function of $(\tilde{y}, \theta)$; hence $\tilde y$ and $\theta$ has a joint normal distribution.

This might also be a standard result but I could not find it proven anywhere. I also think that the justification given is not really accurate, since a bivariate quadratic function takes the form
$$g(x_1, x_2) = Ax_1^2 + Bx_2^2 + Cx_1x_2 + Dx_1 + Ex_2 + F$$
where $A,...,F$ are fixed coefficients and therefore it should be possible that for example $A=B=0$, which surely wouldn't result in a joint normal distribution. But I think what the authors want to say is something like this:
If the integrand is an exponential of a quadratic function on the form
$$g(\tilde{y}, \theta) = -\dfrac{1}{2} \bigg[C_1(x_1 - x_2)^2 + C_2(x_1 - h_1)^2 + C_3(x_2-h_2)^2 + K \bigg] \quad (1.)$$
where $K, h_k$ are constants, $C_k \geq 0$, $C_1 > 0$ and $C_2 > 0$ or $C_3 > 0$, then $(x_1, x_2)$ do have a joint normal distribution (this is just a slight generalization of the example in the book). To see why, we can rewrite $(1.)$ as
$$g(x_1, x_2) = -\dfrac{1}{2} \bigg[Q(x_1, x_2) + L(x_1, x_2) + K_2 \bigg]$$
where
$$
Q(x_1, x_2) = C_1(x_1-x_2)^2 + C_2 x_1^2 + C_3x_2^2\\
L(x_1, x_2) = -2C_2x_1h_1 - 2C_3x_2h_2
$$
and $K_2$ is a constant. Since $C_2 > 0$ or $C_3 > 0$ we have that $Q(x_1, x_2) > 0$ for all $x \neq 0$ so it is a positive definite quadratic form with associated (symmetric) positive definite matrix $A$: $Q(x_1, x_2) = x^TAx$.
The exponent of a multivariate distribution of $x_1, x_2$ can be written as
$$-\dfrac{1}{2}(x^T - \mu^T)P(x-\mu)$$
where $P$ is the precision matrix $\Sigma^{-1}$ and $\mu$ is the mean vector. In order for $(x_1, x_2)$ to have a multivariate normal distribution, we would therefore require that there exists $P, \mu, K_3$ such that
$$Q(x_1, x_2) + L(x_1,x_2) + K_2 =  (x^T - \mu^T)P(x-\mu) + K_3$$
where $K_3$ is a constant.
In other words that there exists $P, \mu$ such that
$$x^TAx + L(x_1, x_2) = x^TPx - 2x^TP\mu \quad (2.)$$
Matching up terms we see that we must have $A = P$, so the question becomes if there exists a mean vector $\mu$ such that
$$L(x_1, x_2) = -2x^TA\mu \iff x^Tu = x^TA\mu \iff u = A\mu,$$
where $u = (C_2h_1, C_3h_2)$. Since $A$ is positive definite it has an inverse so that
$$u = A\mu \iff \mu = A^{-1} u$$
which gives us our required $\mu$. Therefore we see that there does exist $P, \mu$ such that $(2.)$ is satisfied, and therefore $x_1, x_2$ do have a multivariate normal distribution.

We can use the example presented in the book as a check:
In the book example, we have the exponent
$$
-\dfrac{1}{2} \bigg[ C_1(\tilde{y} - \theta)^2 + C_2(\theta - \mu_\Theta)^2 \bigg]
$$
(writing $\mu_\theta$ instead of $\mu_1$ so as to not confuse the notation with the one used above)
where
$$
\begin{align*}
C_1 = \sigma^{-2}\\
C_2 = \tau_1^{-2}
\end{align*}
$$
so that
$$
\begin{align*}
&P = A = \\
&\begin{bmatrix}
C_1 & -C_1\\
-C_1 & C_1^2 + C_2^2
\end{bmatrix} = \\
&\begin{bmatrix}
\sigma^{-2} & -\sigma^{-2}\\
-\sigma^{-2} & \sigma^{-2} + \tau_1^{-2}
\end{bmatrix}
\end{align*}
$$
and from this we get our covariance matrix $\Sigma$ as
$$
\begin{align*}
&\Sigma = P^{-1} = \\
&\begin{bmatrix}
\sigma^{-2} & -\sigma^{-2}\\
-\sigma^{-2} & \sigma^{-2} + \tau_1^{-2}
\end{bmatrix}^{-1} = \\
&\begin{bmatrix}
\sigma^{2} + \tau_1^2 & \tau_1^2\\
\tau_1^2 & \tau_1^{2}
\end{bmatrix}
\end{align*}
$$
so we can see here that $\text{var}(\tilde{Y}|Y) = \sigma^2 + \tau_1^2$ as expected.
The vector $u$ is  $$u = \begin{bmatrix}
0\\
C_2\mu_\Theta
\end{bmatrix} = 
\begin{bmatrix}
0\\
\tau_1^{-2}\mu_\Theta
\end{bmatrix}$$
so that
$$\mu = A^{-1}u = \\
\begin{bmatrix}
\sigma^{2} + \tau_1^2 & \tau_1^2\\
\tau_1^2 & \tau_1^{2}
\end{bmatrix}
\begin{bmatrix}
0\\
\tau_1^{-2}\mu_\Theta
\end{bmatrix} = \begin{bmatrix}
\mu_\Theta\\
\mu_\Theta
\end{bmatrix} $$
so that $E[\tilde{Y}|{Y}] = \mu_\Theta (= \mu_1 \text{ in the book})$ as expected.

I'm sure this is a roundabout way of showing that the authors argument is valid, and there's probably a much more general result (at least I think the argument presented here should be possible to generalize to more than two variables), but I think it is closer to what they had in mind than considering the pdf of the sum of two normal independent random variables.
