Gaussian likelihood + which prior = Gaussian Marginal? Given a Gaussian likelihood for a sample $y$ like $$p(y|\theta) = \mathcal{N}(y;\mu(\theta),\Sigma(\theta))$$ with $\Theta$ being the parameter space and $\mu(\theta)$, $\Sigma(\theta)$ arbitrary parameterisations of the mean vector and the covariance matrix.
Is it possible to specify a prior density $p(\theta)$ and parameterisation of the mean vector $\mu(\theta)$ and the covariance matrix $\Sigma(\theta)$ such that the marginal likelihood $$p(y)=\int_{\theta\in\Theta}N(y;\mu(\theta),\Sigma(\theta))p(\theta)d\theta$$ is a Gaussian likelihood? 
I guess excluding the trivial solution that the covariance is known, that is, $\Sigma(\theta)=\Sigma$, where $\Sigma$ is an arbitrary fixed covariance matrix, this is not possible.
For the special case $\mu(\sigma^2)=\mu$ and $\Sigma(\sigma^2)=\sigma^2$, that is $y$ is one-dimensional, and $p(\sigma^2)=\mathcal{U}(\sigma^2;a,b)$, where $\mathcal{U}(\sigma^2;a,b)$ denotes the uniform density  I can show it: 
\begin{align} 
p(y)&=\int_0^\infty \mathcal{N}(y;\mu,\sigma^2)\mathcal{U}(\sigma^2;a,b)d\sigma^2 \\
&= \frac{1}{b-a} \underbrace{\int_a^b \mathcal{N}(y;\mu,\sigma^2)}_\text{not a Gaussian density}
\end{align}
The accepted answer contains a formal or informal proof or pointers to it.
 A: Your conjecture seems to be true: only a constant variance can lead to
a normal margin. My proof is limited to the case where the
expectation $\boldsymbol{\mu}$ is known, and hence can be assumed to
be zero. For the general case, more sophisticated arguments from
functional analysis seem to be required.
Note that the question is actually about continuous mixture of normals as well
as about Bayes. The statement proved here it that a (continuous) scale
mixture of normals can be normal only for a trivial mixture.
First consider the case of a one-dimensional normal with known mean
$\mu = 0$ and precision parameter $\omega := 1 / \Sigma >0$. Without loss of
generality, we can assume that the parameter $\boldsymbol{\theta}$ is
the precision $\omega$ itself. If the marginal distribution of $y$ is
normal, then $\int \exp\{-y^2 \omega / 2\}\,\omega^{1/2}
p(\omega)\,\text{d}\omega$ is a normal density up to a multiplicative
constant. This density being an even function of $y$ must take the
form $c\exp\{ -y^2 \omega_0 / 2\}$ for some $\omega_0 >0$ and some
constant $c >0$. Since this holds for any $y$ we get with $s := y^2$
$$ \int_0^\infty \exp\{-s \omega \,/ 2\}\,\omega^{1/2}
p(\omega)\text{d}\omega = c \exp\{ -s \omega_0 \,/ 2\} $$ for all $s
\geq 0$, which shows that the finite measure with density function
$\omega \mapsto \omega^{1/2} p(\omega)$ is proportional to the Dirac
mass at $\omega_0$ because these two measures have the same Laplace
transform, up to a multiplicative constant. Thus $\omega$ is almost
surely (a.s.) equal to $\omega_0$.  
This proof extends to the
$d$-dimensional normal with mean zero and precision matrix
$\boldsymbol{\Omega}:=\boldsymbol{\Sigma}^{-1}$. The margin then writes
as $\propto \int \exp\{-\mathbf{y}^\top \boldsymbol{\Omega}\,\mathbf{y} \,/
2\}\,
\left|\boldsymbol{\Omega}\right|^{1/2}p(\boldsymbol{\Omega})\,\text{d}\boldsymbol{\Omega}$
where the integral is on the set $\mathcal{P}$ of positive definite
symmetric $d \times d$ matrices. If this integral is identical to
$c\exp\{ -\mathbf{y}^\top \boldsymbol{\Omega}_0 \mathbf{y} / 2\}$,
then by taking $\mathbf{y}:= \sqrt{s} \,\boldsymbol{u}$ for a scalar
$s \geq 0$ and an arbitrary vector $\mathbf{u}$, we find as above that
$\mathbf{u}^\top \boldsymbol{\Omega}\, \mathbf{u}$ must be a.s. equal
to $\mathbf{u}^\top \boldsymbol{\Omega}_0 \mathbf{u}$, which shows
that $\boldsymbol{\Omega}$ is a.s. equal to $\boldsymbol{\Omega}_0$.
The proof works even if the measure conveniently written as having
density $|\boldsymbol{\Omega}|^{1/2} p(\boldsymbol{\Omega})$
concentrates on a subset of $\mathcal{P}$ with Lebesgue measure zero,
because the Laplace transform argument still applies. So the proof
works for a general parameterisation of the precision (or variance)
matrix.
A: Assume that $\mu$ and $\Sigma$ are a priori independent and that $y$
has a normal margin with mean $\mu_0$ and variance $\Sigma_0$. I
will prove that then the variance $\Sigma$ must be constant, and
the mean $\mu$ must have a normal prior (possibly degenerate).
I will stick to the one-dimensional case for simplicity, using the
characteristic function (c.f.) of $y$, i.e. $\phi_y(t) := 
\mathbb{E}[e^{yit}]$.  We know that $\phi_y(t) = \exp\{\mu_0 it - \Sigma_0 t^2 /2$} and
 a similar formula holds for the distribution of $y$ conditional on $\mu$
and $\Sigma$, which is normal by assumption. So for any real $t$ 
$$
 \mathbb{E}[e^{yit}] = \int \mathbb{E}\left[e^{yit} \, \vert \,\mu,\,\Sigma\right]\,
 p(\mu) p(\Sigma)
 \,\text{d}\mu \text{d} \Sigma =
 \int \exp\left\{ \mu it  - \Sigma t^2/2 \right\} \,p(\mu) p(\Sigma)
 \,\text{d}\mu \text{d}\Sigma, 
$$
and by rearranging the integral, we must have
$$
   \exp\left\{ \mu_0 it - \Sigma_0 t^2 /2 \right\} = 
 \left[\int \exp\left\{ \mu it \right\} p(\mu) 
 \,\text{d}\mu \right]
 \left[\int \exp\left\{ -\Sigma t^2/2\right\} p(\Sigma) 
 \,\text{d}\Sigma \right].
$$
The assumptions needed for such a rearrangement are easily checked.
The first integral at right hand side, say $\phi_1(t)$, is the c.f. of
$\mu$. Note that since $\phi_1(t) e^{-\mu_0 it}$ is found to be real, we see that
the distribution of $\mu$ is symmetric w.r.t. $\mu_0$, and hence that
$\mathbb{E}[\mu] = \mu_0$, as it might have been anticipated.
Now it turns out that the second integral at right hand side, say
$\phi_2(t)$, is also a c.f. To see that, we must check that $\phi_2(0)
= 1$, that $\phi_2$ is continuous at $t=0$ and also that the function
$\phi_2$ is positive definite (p.d.). The first requirement is
obvious, the second is proved by dominated convergence. Now turn to
the p.d. requirement: if the prior distribution written as
$p(\Sigma)\text{d}\Sigma$ is a Dirac mass, then $\phi_2$ is
p.d. because $\phi_2$ is then the c.f. of a normal distribution.  If
the prior is a discrete mixture of Dirac masses, this is true as well
since $\phi_2$ then is the c.f. of a mixture of normals. By a continuity
argument, we see that $\phi_2$ is p.d.
Now let us use the powerful Lévy-Cramér theorem which tells that
both functions $\phi_j$ for $j=1$, $2$ must take the form $\exp\{a_j i t -
b_jt^2 /2 \}$ with $a_j$ real and $b_j \geq 0$. So $\mu$
must be normal (possibly degenerate) with mean $a_1 = \mu_0$.
By simple algebra we then have
$$
    \exp\{ -(\Sigma_0 - b_1) t^2 /2 \} = \int_0^\infty \exp\{ - \Sigma t^2 /2\} p(\Sigma)
    \, \text{d} \Sigma 
$$
which holds for any real $t$. Since any non-negative real writes as $t^2/2$, we see that
the Laplace transform of the prior of $\Sigma$ must be equal to that
of the Dirac mass at $\Sigma_0 - b_1$ and we are done.
A: I have a proposition of proof for you, but you need to check it.
Assume that the marginal likelihood is Gaussian :
$p(y)=\mathcal{N}(y,m,\Gamma)$
then the prior density can be defined by
$p(\theta)=\mathcal{N}(y,\mu(\theta),\Sigma(\theta))^{-1}\mathcal{N}(y,m,\Gamma)f(\theta)$
where $f$ checks $\int_{\theta\in\Theta}f(\theta)d\theta =1$ and $f(\theta)\geq 0$ for $\theta\in\Theta$. ($f(\theta)$ is $p(\theta|y)$).
To be a density, the integral of the prior density $p(\theta)$ on $\Theta$ has to be equal to 1. 
In other words, 
$\int_{\theta\in\Theta}\mathcal{N}(y,\mu(\theta),\Sigma(\theta))^{-1}\mathcal{N}(y,m,\Gamma)f(\theta)d\theta =1$.
It leads to
$\int_{\theta\in\Theta}\mathcal{N}(y,\mu(\theta),\Sigma(\theta))^{-1}\mathcal{N}(y,m,\Gamma)f(\theta)d\theta = \int_{\theta\in\Theta}f(\theta)d\theta$
This equality being true if and only if $\mu(\theta)=m$ and $\Sigma(\theta)=\Gamma$.
