How to solve this integral in moment generating function The moment generating function of generalised Pareto distribution eventually comes down to the following integral (here).
$$
M_X(\theta) = \mathbb Ee^{X\theta} = \int_\mu^\infty e^{\theta x}\frac{1}{\sigma} \left[ 1+\xi\frac{x-\mu}{\sigma} \right]^{-\frac{1}{\xi}-1} dx = e^{\theta\mu} \sum_{j=0}^\infty\frac{(\theta\sigma)^j}{\prod_{k=0}^j(1-k\xi)}, (k\xi < 1),
$$
where $\theta, \xi, \mu \in \mathbb R, \sigma \in \mathbb R^+$. Muraleedharan and Soares (2014) claimed the above integral can be solved by

"adding the integral of each product obtained by multiplying each term of $e^{\theta x}$ with $\frac{1}{\sigma} \left[ 1+\xi\frac{x-\mu}{\sigma} \right]^{-\frac{1}{\xi}-1}$ and substituting $\left[ 1+\xi\frac{x-\mu}{\sigma} \right]^{-\frac{1}{\xi}}=y$ in the integrals of the products".

I do not understand what they are suggesting. Could anyone explain how to do this integral, please? Thank you!
 A: This result just looks wrong to me, like they generalised something beyond its range of validity. What they seem to want done is a power series expansion of the $e^{\theta x}$, followed by some sort of substitution. The problem is that then the swap they sum and the integral, but integrals of the form
$$ \int_a^{\infty} y^n (1+y)^{-a} \, dy $$
don't exist for $n-a \geqslant -1$.
The first obvious thing to do is substitute $\mu + \sigma u= x$, so $du = dx/\sigma$. Then the integral becomes
$$ e^{\theta \mu} \int_0^{\infty} e^{\theta \sigma u} (1+\xi u)^{-1/\xi-1} \, du, \tag{1} $$
which only converges for $\theta \leqslant 0$ (since $\sigma > 0$) and $\xi>0$.
The suggestion appears to now be to look at
$$ \int_0^{\infty} u^j (1+\xi u)^{-1/\xi-1} \, du, $$
and change variables to $ y = (1+\xi u)^{-1/\xi} $. Presumably this is motivated by 
$$ dy = -(1+\xi u)^{-1/\xi-1} \, du $$
At $u=0$, $y=1$, and at $u=\infty$, $y=0$ if $\xi>0$. We also have $u=\frac{1}{\xi}(y^{-\xi}-1)$ so we are left with
$$ \frac{1}{\xi^j} \int_0^1 (y^{-\xi}-1)^{j} \, dy, $$
which does give the
$$ j!\left(\prod_{k=0}^j (1-k\xi)\right)^{-1} $$
answer they have (induction and integration by parts, I suspect). But, this requires $j\xi<1$, or the integrand is too singular at $y=0$ anyway. Hence this method can't possibly derive the whole series legitimately, and the series expansion suggested can't be right.
Indeed, it isn't right: doing the integral (1) exactly gives
$$ e^{\theta (\mu + \sigma/\xi)} \xi^{-1/\xi-1} (-\theta \sigma)^{1/\xi} \Gamma(-1/\xi,-\theta\sigma/\xi), $$
where 
$$\Gamma(a,z)=\int_z^{\infty} t^{a-1} e^{-t} \, dt \tag{2}$$
is the upper incomplete Gamma function. Getting this answer is just a matter of substituting to change (1) into something of the form (2): setting $t=1+\xi u$ will do it.
Suffice to say that this has a power series which contains the given one, but has the extra part
$$ (-\theta )^{1/\xi } \left(-\left(\frac{\sigma }{\xi }\right)^{1/\xi } \Gamma \left(\frac{\xi -1}{\xi }\right)-\theta  \xi ^{-\frac{1}{\xi }-2} \sigma^{\frac{1}{\xi }+1} \Gamma \left(-\frac{1}{\xi }\right)+ \frac{1}{2} \theta ^2 \xi ^{-\frac{1}{\xi }-3} \sigma ^{\frac{1}{\xi }+2} \Gamma \left(-\frac{1}{\xi }\right) + O(\theta^3)\right), $$
which is not taken into account by their answer. Note that what the series in this form tells us is that the derivatives exist in the usual way until the $(-\theta)^{1/\xi}$ appears, and hence only the first few moments (those with $j\xi<1$) exist. This is, of course, exactly what one would expect with a distribution with a density function that behaves as $x^{-1/\xi-1}$ for large $x$.
