Mariusz: I think there are two issues: (1) Despite calling for "numerical" integration, I think a symbolic integration is what is wanted, and (2) The posting is a request as opposed to a question and I can't tell if that is a language issue or a point of view issue.
It appears that what is wanted is a symbolic answer that will give a moment generating function in terms of
$\psi$,
$\phi$, and
$\lambda$.
The density function for the random variable
$Y$ is
f = (1 + (-1 + (1 - Exp[-\[Delta]*y])^-\[Psi])^\[Phi])^-2*(\[Delta] *\[Phi]* \[Psi]*
Exp[-\[Delta]*y]*(1 -Exp[-\[Delta]*y])^(-1 - \[Psi])*(-1 + (1 - Exp[-\[Delta]*y])^-\[Psi])^(-1 + \[Phi]));
$$f=\frac{\delta \psi \phi \exp (-\delta y) (1-\exp (-\delta y))^{-\psi -1} \left((1-\exp (-\delta y))^{-\psi }-1\right)^{\phi -1}}{\left(\left((1-\exp (-\delta y))^{-\psi }-1\right)^{\phi }+1\right)^2}$$
If one sets
$\phi \to 1$, then with there seems to be a symbolic answer with
$\psi>0$ and
$\delta>0$:
mgf = Integrate[Exp[t y] f /. {\[Phi] -> 1}, {y, 0, \[Infinity]},
Assumptions -> {\[Psi] > 0, \[Delta] > 0, t < \[Delta]}]
$$\frac{\psi \Gamma (\psi ) \Gamma \left(1-\frac{t}{\delta }\right)}{\Gamma \left(-\frac{t}{\delta }+\psi +1\right)}$$
I don't know if there is a symbolic (i.e., nice closed form) result
