Fitting data to an Normal Inverse Gaussian distribution

I am using data which I suspect is distributed according to the Normal Inverse Gaussian distribution (NIG). The literature says that the 4 parameters of the NIG distribution can be estimated from the sample mean, variance, kurtosis, and skewness. I would like to test the fit of NIG with the calculated parameters and the data. Now, my function is written as:

NIG2[x_, {α_, β_, μ_, δ_}] := (α/(δ π)) Exp[Sqrt[α^2 - β^2] - (β μ)/δ] (K1[α*Sqrt[1 + ((x - μ)/δ)^2]] Sqrt[1 + ((x - μ)/δ)^2]) Exp[β/δ x]; 

Edit: The K1 function is the Bessel function of the third kind which is:

K1[w_] = Integrate[ Exp[ -w (t + t^(-1)) / 2 ] , {t, 0, Infinity}] (* = ConditionalExpression[2 BesselK[1, w], Re[w] > 0] *) 

The four parameters are

{α, β, μ, δ} 

and their values are (respectively):

{0.635655, - 0.00143495 I, 6.87671, 1.56309} 

I saw Mr Alpha's answer to a similar query where data followed an AR(4) process. I tried using the residualAnalysis module presented, but Mathematica claims that NIG2 is not a valid process. Yet the 4 plots do appear to be quite correct as the residuals seem distributed normally.

How should I test data against a NIG distribution function?