I aim at improving my understanding of Bayesian model averaging in the context of predictions and I am somewhat stuck on how to implement this numerically. Lets say I have two models $\mathcal{M}_1$ and $\mathcal{M}_2$, containing likelihood-functions $\mathcal{L}_i(y_t|\Theta_i)$ and priors for the parameters $\Theta_i$ for the time-series of interest $y_t$. Furthermore I have prior beliefs regarding the models, denoted as $\pi(\mathcal{M}_k)$. After observing $D=(y_1,\ldots,y_T)$ I would like to compute predictions $y_{T+1}$ by evaluating the posterior predictive of the two models. In my understanding this takes the form: $$p(y_{T+1}|D)=\sum_{k=1}^2 \underbrace{p(y_{T+1}|D,\mathcal{M}_k)}_\text{MCMC draws}\underbrace{p(\mathcal{M}_k|D)}_?$$ Given MCMC draws $\{\tilde{y}_{T+1} ^{(i)}\}_{\mathcal{M}_i}$ from the posterior predictive for each model and $p(\mathcal{M}_k|D)$ I would approximate $p(y_{T+1}|D)$ by averaging the draws, weighting with $p(\mathcal{M}_k|D)$ .
According to Bayes Theorem I know that $$p(\mathcal{M}_k|D)\propto p(D|\mathcal{M}_k)\pi(\mathcal{M}_k)$$ where $$p(D|\mathcal{M}_k)=\int p(D|\Theta_k,\mathcal{M}_k)p(\Theta_k|\mathcal{M}_k)d\Theta_k \\ \approx \sum_{i=1}^M p(D|\Theta_k ^{(i)},\mathcal{M}_k)$$ given $M$ MCMC draws from the posterior distribution of the parameters $\Theta_k$.
My question is: Is this correct so far? And if yes, would it be correct to compute $p(D|\mathcal{M}_k)$ based on the approximation first and to weight the draws from the posterior predictive in the next step accordingly?