I am trying to evaluate the following integral marginalized across all possible functions. $$\mathbb{P}(y|X,\theta) = \int \mathbb{P}(y|f)\ \mathbb{P}(f|X,\theta) \ df$$ In G.P. we assume prior to be of zero mean and a covariance matrix i.e. $\mathbb{P}(f|X,\theta) \sim \mathcal{N}(f;0, K)$. The likelihood is modeled as $\mathbb{P}(y|f) \sim \mathcal{N}(y;f,\sigma_n^2I)$ with a mean $f$ and covariance $\sigma_n^2I$.
So we have; $$\begin{align} \mathbb{P}(y|X,\theta) &= \int \mathbb{P}(y|f)\ \mathbb{P}(f|X,\theta)\\ &= \int \mathcal{N}(y;f,\sigma_n^2I)\ \ \mathcal{N}(f;0, K) \tag{1} \end{align}$$
As per http://gaussianprocess.org/gpml/chapters/RW.pdf equation (2.30); the log marginal likelihood is $$\text{log}\ \mathbb{P}(y|X,\theta) = -\frac{1}{2}y^TK_y^{-1}y - \frac{1}{2}\text{log}|K_y| - \frac{n}{2}\text{log}\ 2\pi \tag{2}$$
This means the above integral should give me $\mathcal{N}(y;0,K_y)$ where $K_y = K + \sigma_n^2I$. Furthermore, Appendix (A.7) in the same book gives the product of two Gaussians as: $$\mathcal{N}(x;a,A)\ \mathcal{N}(x;b,B) = z^{-1}\mathcal{N}(x;c,C) \tag{3}$$ where $z^{-1} = \frac{1}{(2\pi)^{n/2}}|A+B|^{-\frac{1}{2}}\text{exp}\big\{-\frac{1}{2}(a-b)^T(A+B)^{-1}(a-b)\big\}$ i.e. $z^{-1} = \frac{1}{(2\pi)^{n/2}}|K_y|^{-\frac{1}{2}}\text{exp}\big\{-\frac{1}{2}f^T(K_y)^{-1}f\big\}$
$c = C(A^{-1}a+B^{-1}b) = (\sigma_n^{-2}I + K^{-1})\sigma_n^{-2}f$
$C = (A^{-1} + B^{-1})^{-1} = \sigma_n^{-2}I + K^{-1}$
With this in mind, if I plug (3) in (1) and carry out the integration I should simply get $z^{-1}$ with $f$ replaced by $y$ (not quite sure if I can do this) and since $\mathcal{N}(c,C)$ is a pdf it should give me a value of 1. If I do this I end up with (2) after taking log. But I am not 100% confident on my approach.
Any help in solving the integral would be appreciated.