Cumulative Distribution Function of the skewed generalized error distribution for qq-plot

I´m trying to calculate the cumulative distribution function of the skewed generalized error distribution with the probability density function from Theodossiou(http://www.mfsociety.org/modules/modDashboard/uploadFiles/journals/MJ~0~p1a4fjq38m1k2p45t6481fob7rp4.pdf):

And in R it looks like this:

psi <- -0.09547862
m <- 0.1811856
g <- -0.1288893
d <- 0.8029088

c <- (2/(1+exp(-g)))-1
p <- exp(psi)

y <- function(x) ((d**(1-(1/d)))/(2*p))*gamma(1/d)**(-1)*exp(-(1/d)*((abs(x-m)**d)/((1+sign(x-m)*c)**(d)*p**(d))))

I do this is to fit the skewed generalized error distribution to my data and assess the distributions fit to my data by creating a qq-plot. So now I need to calculate the cumulative distribution function and then the inverse cdf. For the inverse cdf I plan to use the inversion-function from the GofKernel-Package. But for this I need the cdf. Is there anyway to calculate that with numerical integration in R?

To get a cumulative function via integration you can pass the x-values to a function that integrates from a suitable extreme low value to an upper limit that is x

 # First look at the density function
 plot( y(x) ~ x )


 cum <- sapply(x, function(x) integrate(y,-10, x)$value )
 plot( cum ~ x)
 
 # So the inverse is just `x`  as a function of `cum`
 plot( x ~ cum)