On the computation of the cumulative distribution function of the Normal Inverse Gaussian distribution

In this paper, we obtain various series and asymptotic expansions involving the modified Bessel function of the second kind for the normal inverse Gaussian cumulative distribution function. The new expansions accelerate computations, complementing the numerical integration methods implemented in statistical software packages. We also provide a detailed description of the algorithm and its corresponding implementation in C++. The performance and accuracy of the algorithm are extensively tested and benchmarked with open-source implementations, offering superior accuracy and speed-ups of a factor from 5 to 60.