I am trying to solve the integral
$$I(y)=\int_{\mathbb R}f(x,y)g(x)\,dx,$$
where $f(x,y)$ is the bivariate normal density with known mean $(\mu_1,\mu_2)$ and covariance matrix $\Sigma$ , and $g(x)$ is a normal density with mean $\mu$ and variance $\sigma^2$. I am stuck on this and would appreciate some help (I guess there is some sort of trick to do this without too much pain).
In case someone is interested in the motivation: this integral appears in a measurement error model where one observes $(x,y)$ which is a noisy observation of the quantities of interest $(\mu_1,\mu_2)$.
A common approach is to express the integrand in the form $\exp(ax^2+bx+c)$ and then complete the square in the argument of the exponential. Then, what remains can be expressed in the form of an integral of a normal density over $\mathbb R$ which has value $1$, and of course anything involving $y$ and various constants is what is left. In short, $I(y)$ should work out to be of the form $\alpha\exp(-\beta (y-\mu_Y)^2)$. For example, $$\begin{align} \exp(-(x^2 + y^2 -xy)) &= \exp(-(x^2 - 2x(y/2) + y^2/4 + 3y^2/4))\\ &= \exp(-3y^2/4)\exp(-(x-y/2)^2) \end{align}$$ where $\exp(-3y^2/4)$ and sundry constants can be pulled outside the integral, leaving the integral of a normal density with mean $-y/2$ and variance $\ldots$, which integrates to $1$. So the whole thing simplifies to an exercise in high-school algebra.
Since the OP requested "a more elegant trick", here goes. $f_{X,Y}(x,y) = f_{X\mid Y}(x\mid y)f_Y(y)$ where $f_{X\mid Y}(x\mid y)$ and $f_Y(y)$ are $N(\mu_1+\rho\frac{\sigma_1}{\sigma_2}(y-\mu_2),(1-\rho^2)\sigma_1^2)$ and $N(\mu_2,\sigma_2^2)$ random variables where $\rho$ is the correlation coefficient of $X$ and $Y$.
Fix $y$ to be some constant, and consider a normal random variable $\hat{X}$ with mean $\mu_1+\rho\frac{\sigma_1}{\sigma_2}(y-\mu_2)$ and variance $(1-\rho^2)\sigma_1^2$. Let $Z$ be a $N(\mu,\sigma^2)$ random variable that is independent of $\hat{X}$ and note that $f_Z(x) = g(x)$. Then $W = \hat{X}-Z$ is also a normal random variable with mean $\mu_W = \mu_1+\rho\frac{\sigma_1}{\sigma_2}(y-\mu_2) - \mu$ and variance $\sigma_W^2 = (1-\rho^2)\sigma_1^2 + \sigma^2$ so that $$f_W(0) = \frac{1}{\sigma_W\sqrt{2\pi}}\exp\left(-\frac{\mu_W^2}{2\sigma_W^2}\right).$$ But, $$\begin{align} f_W(0) &= \int_{-\infty}^\infty f_{\hat{X}}(x)f_Z(x)\,\mathrm dx\\ &= \int_{-\infty}^\infty \frac{f_{X,Y}(x,y)}{f_Y(y)} g(x)\,\mathrm dx\\ \Rightarrow \qquad f_W(0)f_Y(y)&= \int_{-\infty}^\infty f_{X,Y}(x,y) g(x)\,\mathrm dx\\ \Rightarrow \qquad I(y) &= f_W(0)f_Y(y)\\ &= \frac{1}{\sigma_W\sqrt{2\pi}}\exp\left(-\frac{\mu_W^2}{2\sigma_W^2}\right) \frac{1}{\sigma_2\sqrt{2\pi}}\exp\left(-\frac{(y-\mu_2)^2}{2\sigma_2^2}\right). \end{align}$$ Now, plug-and-chug for $\mu_W$ and $\sigma_W$ to get the final result.
Is that any easier or any more elegant?