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)

In general, if you want to estimate the cumulative distribution function, use the function ecdf as follows:

x <- seq(-10,10,0.1)
Fn <- ecdf(y(x))
plot(Fn)

If you want to visualize how two data sets are similar, use qqplot as follows:

y1 <- y(x)            # from your function
y2 <- rnorm(100)      # some generic data
qqplot(y1, y2)        # if the two data sets are from the same 
                      # distribution, you should see a straight line