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!
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$.