Creating values from a normal inverse gaussian (NIG) distribution

I have a vector of empirical observations ‘a’ containing about 16000 values.

I’m trying to test, whether the normal inverse gaussian (NIG) fits the empirical data or not. First, I estimate the distribution parameters on the data:

SP500NIG<-fit.NIGuv(a, silent = TRUE)

Then I try to generate values from the NIG distribution with the parameters above. I do the following:

library(fBasics)

SP_NIG_sim<-rnig(length(a),mu=SP500NIG@mu,delta=SP500NIG@sigma,alpha = SP500NIG@alpha.bar,beta=SP500NIG@gamma)

And get the corresponding theoretical histogram, which doesn’t look like the empirical one (left - empirical one, right - theoretical one):

Where is the mistake in my code? How to generate values from the NIG distribution with the estimated parameters, which will more or less coincide with the empirical data?

Thank you for your help.