Is there an elegant (probabilistic) way to compute $E(X^n e^{-\lambda X})$, where $n \in \mathbb{N}$, $\lambda > 0$ and $X$ is a random variable with normal distribution $N(\mu,\sigma^2)$?
Alternatively, my question is simply how to compute the integral $$ \frac{1}{\sqrt{2\pi\sigma^2}} \int_{-\infty}^{\infty} x^n e^{-\lambda x} e^{-\frac{(x-\mu)^2}{2\sigma^2}} dx. $$
On
HINT You have $$ \begin{split} \lambda x + \frac{(x - \mu)^2}{2\sigma^2} &= \frac{x^2 + \mu^2 + 2x(\sigma^2\lambda-\mu)}{2\sigma^2} \\ &= \frac{x^2 + 2x(\sigma^2\lambda-\mu) +(\sigma^2\lambda-\mu)^2} {2\sigma^2} + \frac{\mu^2 - (\sigma^2\lambda-\mu)^2}{2\sigma^2}\\ &= \frac{(x + \sigma^2\lambda-\mu)^2} {2\sigma^2} + \frac{\mu^2 - (\sigma^2\lambda-\mu)^2}{2\sigma^2}\\ \end{split} $$ and denote the last term on the RHS by $K$ and you have $$ \int_\mathbb{R} x^n e^{\lambda x} e^{\frac{(x - \mu)^2}{2\sigma^2}} dx = e^K\int_\mathbb{R} x^n \exp\left(\frac{(x + (\sigma^2\lambda-\mu))^2}{2\sigma^2}\right) dx = e^K \mathbb{E}[Y] $$ where $Y$ is also normal with same std deviation but shifted mean.
On
The exponential generating function of the sequence $A_n = \mathbb E[X^n e^{-\lambda X}]$ is $$ g(z) = \sum_{n=0}^\infty \frac{A_n}{n!} z^n = \mathbb E[e^{(z-\lambda)X}] = \exp\left(\frac{\sigma^2 (z-\lambda)^2}{2} + 2 \mu (z - \lambda) \right)$$ and $A_n$ can be obtained from the Maclaurin series of that.
On
First, let $$ f(x) = \dfrac{1}{\sigma\sqrt{2\pi}}e^{-\frac{1}{2}\left(\frac{x - \mu}{\sigma}\right)^2}, $$ and recall that the moment generating function of the normal distribution with mean $\mu$ and variance $\sigma^2$ is given by $$ M_X(t) = E\left(e^{tX}\right) = \exp\left(\mu t + \frac{1}{2}\sigma^2 t^2\right). $$ Second, notice that $$ \dfrac{\partial}{\partial \lambda} e^{-\lambda x} = -x e^{-\lambda x}, $$ and $$ \dfrac{\partial^2}{\partial \lambda^2} e^{-\lambda x} = x^2 e^{-\lambda x}; $$ accordingly, by mathematical indication, we get $$ \dfrac{\partial^n}{\partial \lambda^n} e^{-\lambda x} = (-1)^n x^n e^{-\lambda x}. $$ Hence, by using Leibniz integral rule, the given expected value can be rewritten as $$ E\left(X^n e^{-\lambda X}\right) = (-1)^{-n} \int_\mathbb{R} \dfrac{\partial^n}{\partial \lambda^n} e^{-\lambda x} f(x) dx = (-1)^{-n} \dfrac{d^n}{d \lambda^n} \int_\mathbb{R} e^{-\lambda x} f(x) dx = (-1)^{-n} \dfrac{d^n}{d \lambda^n} M_X(-\lambda) = (-1)^{-n} \dfrac{d^n}{d \lambda^n} \exp\left(\frac{1}{2}\sigma^2 \lambda^2 - \mu \lambda\right). $$
Finally, depending on the value of $n$ and $\lambda$, you can determine the desired expectation. For instance, if $n = 1$, then $$ E\left(X e^{-\lambda X}\right) = \left(\mu - \lambda \sigma^2\right) M_X(-\lambda). $$
$$\lambda x + \frac{(x-\mu)^2}{2\sigma^2} = \frac{x^2 - 2(\mu - \sigma^2 \lambda)x + \mu^2}{2\sigma^2} = \frac{(x-(\mu-\sigma^2 \lambda))^2}{2\sigma^2} + \frac{2 \lambda \mu -\sigma^2 \lambda^2}{2}$$
Ignoring constant factors, your integral becomes $$\int x^n e^{-\frac{(x-(\mu-\sigma^2 \lambda))^2}{2\sigma^2}} \mathop{dx},$$ which can be computed from the [non-central] moments of a $N(\mu-\sigma^2 \lambda, \sigma^2)$ distribution.