Is it always better to average out parameter uncertainty? Setup
If we have a data set $y_1, \ldots, y_t := y_{1:t}$, and we're trying to predict $y_{t+1}$, a Bayesian would try to use the posterior predictive distrbution
$$
p(y_{t+1} \mid y_{1:t}) = \int p(y_{t+1} | y_{1:t}, \theta)p(\theta | y_{1:t})d\theta .
$$
as opposed to plugging in some estimator of the parameter into the conditional likelihood:
$$
p(y_{t+1} \mid y_{1:t}, \hat{\theta}).
$$
When asked why the former is superior to the latter, people often reply "because it takes into account parameter uncertainty."

Question
Now, assuming the first is unavailable, and you only care about prediction score, is it always better to average out parameter uncertainty? For instance, what if we average over with the prior distribution:
$$
\check{p}(y_{t+1} \mid y_{1:t}) := \int p(y_{t+1} \mid y_{1:t}, \theta)p(\theta) d\theta?
$$
Letting $f(y_{t+1})$ be the "true" prediction density; can we say anything about information theoretic quantities like
$$
\text{KL}\left[ f(y_{t+1}) || \check{p}(y_{t+1} \mid y_{1:t})\right] \overset{\text{?}}{\le} 
\text{KL}\left[ f(y_{t+1})  || p(y_{t+1} \mid y_{1:t}, \hat{\theta})\right]?
$$
I doubt the above is true for any $\theta$. I can show it's true on average:
$$
\text{KL}\left[ f(y_{t+1}) || \check{p}(y_{t+1} \mid y_{1:t})\right] \le 
E_{\theta} \left\{ \text{KL}\left[ f(y_{t+1})  || p(y_{t+1} \mid y_{1:t}, \theta )\right] \right\}.
$$
Proof: by Jensen's $\log\check{p}(y_{t+1} \mid y_{1:t}) \ge E_{\theta}\left[\log p(y_{t+1} \mid y_{1:t}, \theta) \right]$. Take $E_{f(y_{t+1})}(\cdot)$ on both sides, and then switch the order of integration.
 A: Okay, there are several points in your question, I'll try to address them one-by-one.
First: should you average $\theta$ out if you only care about predicting $y_{t+1}$?
Let's assume we have some joint data-generating process $p(y_{t+1}, y_{1:t}, |\theta)$ and a prior $p(\theta)$. Since you "only care about prediction score" and wanna be Bayesian, you should employ the Bayesian Decision Theory:

*

*You'll have a certain "rule" $a(\cdot)$ that takes $y_{1:t}$ and produces a best guess of what $y_{t+1}$ might be: $$ \hat{y}_{t+1} = a(y_{1:t}) $$

*You need to evaluate how good a particular prediction is. To do so you introduce a loss function $L(\hat{y}, y)$

*There's a lot of uncertainty ($y_{t+1}$ might not be uniquely determined by $y_{1:t}$, and all $y$s depend on the unknown $\theta$), so we average our loss: $$ \mathcal{L}(a) = \mathop{\mathbb{E}}_{p(y_{1:t+1}, \theta)} \; L(a(y_{1:t}), y_{t+1}) = \int L(a(y_{1:t}), y_{t+1}) p(y_{1:t+1}, \theta) dy_{1:t+1} d\theta $$

*Next, you seek to find an optimal decision rule $a(\cdot)$ that minimizes this expected loss.

Here's an illustrative special case. Consider $L^{(2)}(\hat{y}, y) = (\hat{y} - y)^2$. Then
$$ \mathcal{L}^{(2)}(a) = \mathop{\mathbb{E}}_{p(y_{1:t+1}, \theta)} \; (a(y_{1:t}) - y_{t+1})^2 = \mathop{\mathbb{E}}_{p(y_{1:t})} \mathop{\mathbb{E}}_{p(y_{t+1}, \theta \mid y_{1:t})} (a(y_{1:t}) - y_{t+1})^2 $$
It's known that expected L2 loss is minimized by the average of that distribution: $$ a = \mathop{\arg\min}_a \mathbb{E}_{p(x)} (a-x)^2 \Rightarrow a = \mathbb{E}_{p(x)} x $$
Hence, the optimal decision rule for our $L^{(2)}$ loss is:
$$
\begin{align*}
a(y_{1:t})
&= \int\int y_{t+1} p(y_{t+1}, \theta \mid y_{1:t}) dy_{t+1} d\theta \\
&= \int\int y_{t+1} p(\theta \mid y_{t+1}, y_{1:t}) p(y_{t+1} \mid y_{1:t}) dy_{t+1} d\theta \\
&= \int y_{t+1} \left(\int p(\theta \mid y_{t+1}, y_{1:t}) d\theta\right) p(y_{t+1} \mid y_{1:t}) dy_{t+1} \\
&= \int y_{t+1} p(y_{t+1} \mid y_{1:t}) dy_{t+1}
\end{align*}
$$
We've just integrated out $\theta$. Indeed, you can see that since $\theta$ is unknown a-priory we can't use it in the decision rule $a(\cdot)$ and neither do we care about its particular value. This means we could have integrated $\theta$ out long time ago, right in the general expected loss: $$ \mathcal{L}(a) = \mathop{\mathbb{E}}_{p(y_{1:t+1})} \; L(a(y_{1:t}), y_{t+1}) = \int L(a(y_{1:t}), y_{t+1}) p(y_{1:t+1}) dy_{1:t+1} $$
So all this lets us conclude that the optimal decision rules $a(\cdot)$ are some statistics of the posterior-predictive distribution $p(y_{t+1} \mid y_{1:t})$. For squared-difference loss its posterior mean, for absolute-error loss it's posterior median, etc. This shows that, yes, you should average $\theta$ out if you're operating in the Bayesian paradigm. Ultimately, it minimizes average loss averaged over all possible environments ($\theta$s) and outcomes in these environments ($y$s). If your environment is actually fixed once and for all and the prior is wide (meaning, it spreads over many different values of $\theta$), performance of Bayesian procedure will be hindered by this mismatch.
Huh, what a mouthful! Now on, to the second question.

Second: Posteriors are hard! Can I just integrate over the prior $p(\theta)$ instead of the posterior $p(\theta \mid y_{1:t})$?
In the first part we've seen that using the posterior-predictive distribution is the best thing to do (under certain assumptions). Opting for prior-predictive distribution $\check{p}(y_{t+1} \mid y_{1:t})$ is thus suboptimal. But how much?
Well, the answer depends on a particular model $p(y_{1:t+1}, \theta)$. Consider two extreme cases:

*

*First, generate $y_1$ from some simple distribution, say, standard Gaussian $p(y_1) = \mathcal{N}(y_1 \mid 0, 1)$. Next, put $y_{k+1} = y_k$ for $k = 1 \dots t$, and draw $\theta$ from $\mathcal{N}(y_1, 100)$ – a Gaussian with mean $y_1$ and huge variance of 100. You can see that in this model you can predict $y_{t+1}$ quite easily just by copying $y_t$ and the posterior is not needed at all! This is because all the information that's required to construct $y_{t+1}$ is already there in $y_1$ and $\theta$ does not bring anything to the table ($p(y_{t+1} \mid y_{1:t}, \theta)$ does not actually depend on $\theta$ except for $y_1$).

*If we instead generate $\theta$ from some standard Gaussian: $p(\theta) = \mathcal{N}(\theta | 0, 1)$ and then put $y_k = k \theta$ for all $k$, we'll have a different setting. Now, $p(y_{t+1} \mid y_{1:t}, \theta)$ does not use $y_{1:t}$ and all the required information is contained in the $\theta$. Now, it's crucial to perform accurate inference on $\theta\mid y_{1:t}$ – by doing so we essentially "extract" information about $\theta$ from the observed $y_{1:t}$.

These two examples show that there are two information pathways in your model: from shared parameter $\theta$ to each $y_t$ and from the previous element $y_t$ to the next one $y_{t+1}$. If the shared parameter influences each $y_t$ in a non-trivial manner, using a prior-predictive distribution will likely result in a poor approximation. If most of the information is contained in the previous terms, and the shared parameter contributes little, then you'll be fine (although perhaps there's not much sense in introducing the unobserved latent parameter in the first place!).
Finally, you might attempt quantifying how much the prior-predictive distribution worse by calculating the expected loss for a decision rule given by, say, prior-predictive mean: $$ \check{a}(y_{1:t}) = \int y_{t+1} p(y_{t+1} \mid y_{1:t}, \theta) p(\theta) d\theta $$

Third: can a prior-predictive distribution $\check{p}(y_{t+1} | y_{1:t})$ be closer (in terms of KL divergence) to the posterior-predictive $p(y_{t+1} | y_{1:t})$ than a point estimate $p(y_{t+1} | y_{1:t}, \hat{\theta})$ for the best possible $\hat\theta$?
The second example (and its particular model!) from the previous section shows us that if we don't perform posterior inference (that is, don't extract information about $\theta$ from $y_{1:t}$), we won't be able to produce a coherent $y_{t+1}$. Therefore, in this case prior-predictive would be a poor approximation to the posterior-predictive, whereas the point-estimate should be a much better one (especially if the true posterior is sharply concentrated in a small area). On the other hard, in the case of the first example, we'll have both prior predictive and the point estimate would be both good at approximating (perhaps even exactly) the posterior-predictive.
Is the prior-predictive always worse then? I have one more example: consider a multivariate standard Gaussian distribution $\theta \sim \mathcal{N}(0, I)$ of as many dimensions as needed. Then we generate $y$s as $y_k = \theta_k$. That is, we essentially observe $\theta$'s components one-by-one. Obviously, since all components of $\theta$ are independent, the posterior-predictive will be equal to the prior equal to the standard univariate Gaussian. However, the distribution $p(y_{t+1}|y_{1:t}, \theta)$ in this model is degenerate (since $y_{k}$ is a deterministic function of $\theta$) and is given by a delta function $\delta(y_{t+1} - \theta_{t+1})$. Therefore any point-estimate $p(y_{t+1}|y_{1:t}, \hat\theta)$ will also be a degenerate distribution and would be an infinitely worse approximation to the posterior-predictive!
So neither relation holds in a general case.
