integrate product of two gaussian pdfs where one marginal is not gaussian

Question

I'm seeing some help as of how to integrate the product of two Gaussian PDFs analytically, $g(y_n | x_n)$ and $f(x_n | x_{n-1})$ as described by Whiteley (http://www.maths.bris.ac.uk/~manpw/apf_chapter.pdf). \begin{equation} \int\limits_{-\infty}^{\infty}g(y_n | x_n)f(x_n | x_{n-1})dx_n \end{equation} It's well known that the product of two Gaussians $f(x)g(x)$ is a new unscaled Gaussian, however for my case the conditional probability $g(y_n | x_n)$ is Gaussian but the marginal $g(y_n)$ itself may not be Gaussian. That is, we have two time series defined by means $bx_n, ax_{n-1}$ and variances $\sigma_f^2, \sigma_g^2$ respectively (where $\epsilon$ represents white noise): \begin{align} x_n &= ax_{n-1} + \sigma_f\epsilon_f \\ y_n &= bx_n + \sigma_g\epsilon_g \end{align} and I want to integrate the product of the two corresponding PDFs \begin{align} g(y_n | x_n) &=\frac{1}{\sqrt{2\pi}\sigma_g}exp\big(-\frac{(y_n-bx_n)^2}{2\sigma_g^2} \big) \\ f(x_n | x_{n-1} ) &= \frac{1}{\sqrt{2\pi}\sigma_f}\big(-\frac{(x_n-ax_{n-1})^2}{2\sigma_f^2} \big) \end{align} The derivation for the general case is outlined step by step (http://www.tina-vision.net/docs/memos/2003-003.pdf), however for my case the details get a bit tricker. Taking the product of f and g: \begin{equation} f(x)g(x) = \frac{1}{2\pi\sigma_f\sigma_g}exp\Big(-\frac{(y_n-bx_n)^2}{2\sigma_f^2} + \frac{(x_n-ax_{n-1})^2}{2\sigma_g^2} \Big) \end{equation} Expand the two quadratic terms in the exponents and collect terms: \begin{align} Let \beta :&= \Big(-\frac{(y_n-bx_n)^2}{2\sigma_f^2} + \frac{(x_n-ax_{n-1})^2}{2\sigma_g^2} \Big) \\ &= \frac{y_n^2\sigma_g^2 - 2by_nx_n\sigma_g^2 + b^2x_n\sigma_g^2 + x_n^2\sigma_f^2 - 2ax_nx_{n-1}\sigma_f^2 + a^2x_{n-1}^2\sigma_f^2}{2\sigma_f^2\sigma_g^2} \\ &= \frac{y_n^2\sigma_g^2 + x_n^2\sigma_f^2 - 2x_n(by_n\sigma_g^2 - ax_{n-1}\sigma_f^2) + b^2x_n\sigma_g^2 + a^2x_{n-1}^2\sigma_f^2}{2\sigma_f^2\sigma_g^2} \end{align} Is there is a neat way to make this work out nicely to get a closed form solution for the mean, variance and scaling factor of the new Gaussian? Any help would be appreciated!