I would like to know if there's a way to compute or approximate the following expectation:
$$\mathbb{E}[(c+e^X)^{-n}]$$
where $X=N(0,\sigma^2)$ and $n,c>0$ (you can also assume that $n$ is a natural number if necessary). After googling a bit I found the following articles: 1, 2 which talk about negative moments of positive random variables but applying the methods yields a divergent series.
Since you say you're interested in computing the expectation, here's a way to calculate the moments numerically.
Rearranging a little and defining $Y = X - \log c$, the expectation you're after, which I'll call $I$, is
$$I = c^{-n} \cdot \mathbb{E} \left[ \frac{ 1 }{ \left( 1 + e^{Y} \right)^{n}} \right].$$
Massaging the term inside the integral still further, this can be rewritten
$$ I = c^{-n} \cdot \mathbb{E} \frac{ e^{-nY/2} }{ \left( e^{Y/2} + e^{-Y/2} \right)^{n} } = c^{-n} \cdot \mathbb{E} \left[ e^{-nY/2} \frac{ 1 }{ \left( 2 \cosh (Y/2) \right)^{n} } \right] \ .$$
The point of doing this is that $1/(2\cosh Y/2)^{n}$ has a nice Fourier transform. Specifically, $$\frac{1}{(2 \cosh (Y/2))^{n}} = \int_{-\infty}^{\infty} e^{iYz} \mathscr{F}_{n}(z) \, dz \ , $$ where $$ \mathscr{F}_{n}(z) = \frac{1}{2 \pi} \frac{ \Gamma(n/2+iz)\Gamma(n/2-iz)}{ \Gamma(n) } $$ are rapidly decaying, even functions of $z$. (For integer $n$, they can be expressed more simply: for example, $\mathscr{F}_{1}(z) = \frac{1}{2} \operatorname{sech} \pi z$.)
So $$I = c^{-n} \mathbb{E} \int_{-\infty}^{\infty} e^{-n Y/2} e^{iYz} \mathscr{F}_{n}(z) \, dz \ . $$
Writing out the characteristic function of the Normal random variable $Y$ explicitly, and rearranging a bit, this becomes
$$I = c^{-n/2} \int_{-\infty}^{\infty} \cos \left[ \left( \log c + \frac{1}{2} \sigma^{2} n \right) z \right] e^{ \frac{1}{2} \sigma^{2} \left( \frac{n^{2}}{4} - z^{2} \right) } \mathscr{F}_{n}(z) \, dz \ , $$
which is easy to evaluate numerically. Here's what the moments look like for $\sigma=1$ and $c$ equal to 0.7, 1, and 1.3:
This approach is also useful if $X$ and $Y$ are not Normal, so long as you know their characteristic functions.
As $c \to \infty$ for fixed $n$ and $\sigma$, it looks to me like $ {\mathbb E}\left[(c+e^X)^{-n}\right]$ has the asymptotic series $$c^{-n} + \sum_{j=1}^\infty (-1)^j {{n+j-1} \choose j} c^{-n-j} e^{j^2 \sigma^2/2}$$
Yes, like most asymptotic series, this is divergent, but partial sums might still be useful approximations.
Maple 2015 produces a simple answer by
with(Statistics): X := RandomVariable(Normal(0, sigma)):
Mean(c+exp(X))^(-n) assuming c>0, n::posint;
$$ \left( c+{{\rm e}^{1/2\,{\sigma}^{2}}} \right) ^{-n} $$
