Integral of a product between two normal distributions and a monomial The integral of the product of two normal distribution densities can be exactly solved, as shown here for example.
I'm interested in compute the following integral (for a generic $n \in \mathbb{N}$):
$$
I_n = \int_{-\infty}^{\infty} \frac{1}{\sqrt{2 \pi} \sigma} e^{-\frac{x^2}{2 \sigma^2}} \frac{1}{\sqrt{2 \pi} \rho} e^{-\frac{(x-\mu)^2}{2 \rho^2}} x^n\mathrm{d}x
$$
Is it possible to extend this result and solve $I_n$ for a generic $n \in \mathbb{N}$?
 A: Yes, there are closed form expressions for
$$I_n = \int_{-\infty}^{\infty} \frac{1}{\sqrt{2 \pi} \sigma} e^{-\frac{x^2}{2 \sigma^2}} \frac{1}{\sqrt{2 \pi} \rho} e^{-\frac{(x-\mu)^2}{2 \rho^2}} x^n\mathrm{d}x,\;\;n\in\mathbb{N}.$$
For $\mu=0$ this is simplest, one has
$$I_n=(2\pi)^{-1/2}(n-1)\text{!!} \left(\rho^2 \sigma^2\right)^{n/2} \left(\rho^2+\sigma^2\right)^{-n/2-1/2},\;\;\text{for}\;\;n\;\;\text{even},$$
$I_n=0$ for $n$ odd.
For $\mu\neq 0$ the expressions are a bit more lengthy,
$$I_n=Z_nP_n\left(\frac{\mu \sigma}{\rho \sqrt{\rho^2+\sigma^2}}\right),$$
$$Z_n=(2\pi)^{-1/2}e^{-\frac{\mu^2}{2 \left(\rho^2+\sigma^2\right)}} (\rho \sigma)^n \left(\rho^2+\sigma^2\right)^{-n/2-1/2}.$$
The function $P_n(x)$ is a polynomial of degree $n$, the first eight are given by
$$\left\{x,x^2+1,x^3+3 x,x^4+6 x^2+3,x^5+10 x^3+15 x,x^6+15 x^4+45 x^2+15,x^7+21 x^5+105 x^3+105 x,x^8+28 x^6+210 x^4+420 x^2+105\right\}.$$
This is a polynomial described at https://oeis.org/A100861, given by
$$P_n(x)=x^n+\sum _{k=1}^{n/2} \frac{n! x^{n-2 k}}{2^k k! (n-2 k)!}.$$
A: $\def\m{\mu}
\def\p{\pi}
\def\s{\sigma}
\def\f{\varphi}
\def\r{\rho}
\def\mm{M}
\def\ss{S}$Let
\begin{align*}
\f(x;\m,\s) &= \frac{1}{\sqrt{2\p}\s} e^{-(x-\m)^2/(2\s^2)}
\end{align*}
and
\begin{align*}
\f_m(x) &= \prod_{i=1}^m \f(x;\m_i,\s_i) \\
&= \frac{1}{(2\p)^{m/2}\prod_{i=1}^m \s_i}
\exp\left(-\sum_{i=1}^m \frac{(x-\m_i)^2}{2\s_i^2}\right).
\end{align*}
By completing the square one finds
\begin{align*}
\f_m(x) &= A(m)\f(x;\mm,\ss),
\end{align*}
where
\begin{align*}
\frac{1}{\ss^2} &= \sum_{i=1}^m \frac{1}{\s_i^2} \\
\mm &= \sum_{i=1}^m \frac{\m_i}{\s_i^2} \\ 
A(m) &= \frac{\ss}{(2\p)^{(m-1)/2} \prod_{i=1}^m \s_i}
\exp
\left[
\frac12\left(
\frac{\mm^2}{\ss^2} - \sum_{i=1}^m \frac{\m_i^2}{\s_i^2} 
\right)
\right].
\end{align*}
That is, such a product is itself proportional to a normal distribution.
Thus,
\begin{align*}
\int_{-\infty}^\infty x^n \f_m(x) \,dx 
&= A(m) \int_{-\infty}^\infty x^n \f(x;\mm,\ss)\,dx \\ 
&= A(m) \int_{-\infty}^\infty ((x-\mm)+\mm)^n \f(x;\mm,\ss)\,dx \\ 
&= A(m) \int_{-\infty}^\infty 
\sum_{k=0}^n \binom{n}{k} (x-\mm)^k \mm^{n-k} 
\f(x;\mm,\ss)\,dx \\ 
&= A(m) 
\sum_{k=0}^n \binom{n}{k} 
\mm^{n-k} 
\int_{-\infty}^\infty 
(x-\mm)^k 
\f(x;\mm,\ss)\,dx \\ 
&= A(m)\sum_{k=0\atop k{\textrm{ even}}}^n
\binom{n}{k} \mm^{n-k} \ss^k (k-1)!!,
\end{align*}
and so
\begin{align*}
\int_{-\infty}^\infty x^n \f_m(x) \,dx
&= 
\frac{\ss \mm^n}{\prod_{i=1}^m \s_i} 
\frac{1}{(2\p)^{(m-1)/2}}
\exp
\left[
\frac
12\left(
\frac{\mm^2}{\ss^2} - \sum_{i=1}^m \frac{\m_i^2}{\s_i^2} 
\right) 
\right] \\
& \quad \times \sum_{k=0}^{\lfloor n/2\rfloor} 
\binom{n}{2k} (2k-1)!! 
\left(\frac{\ss}{\mm}\right)^{2k} 
\end{align*}
For the original problem,
\begin{align*}
m &= 2 \\ 
(\m_1,\s_1) &= (0,\s) \\
(\m_2,\s_2) &= (\m,\r)
\end{align*}
so
\begin{align*}
\mm &= \frac{\m\s^2}{\s^2+\r^2} \\ 
\ss^2 &= \frac{\s^2\r^2}{\s^2+\r^2}. 
\end{align*}
With a little work, we find a formula agreeing with that of @CarloBeenakker.
