Computing the distribution function of the NIG for small delta

Accurate computation of the cumulative distribution function of the normal inverse Gaussian distribution becomes particularly challenging when both $\delta$ and $|x-\mu|$ are small. In Navas-Palencia (2026), we developed a robust algorithm combining series expansions, asymptotic expansions, and numerical integration. Although the algorithm remains applicable as $\delta\to0$, it must resort to numerical integration in this critical regime, which can reduce performance. For the symmetric case $\beta=0$, we also developed an exponentially improved asymptotic expansion to accelerate the series as $|x-\mu|\to0$. While effective in this regime, that expansion is considerably more involved to implement in the C++ library and has so far only been tested in arbitrary precision using mpmath.

Here, we follow a different approach. We directly expand the factor $e^{-\delta^2/(2t)}$, which leads to a series whose coefficients can be computed efficiently by recurrence. In this post, we focus on the special case $\beta=0$ for brevity and compare the resulting series with the convergence-acceleration method described in Navas-Palencia (2026), Section 3.1.2. Extending the approach to $\beta\neq0$ is straightforward; only the first coefficient, $I_0$, requires additional analysis. This extension is part of ongoing work on the generalized hyperbolic distribution. Finally, the notation in this post follows the, apparently, more common ordering $(\alpha, \beta, \delta, \mu)$, whereas the paper uses $(\alpha, \beta, \mu, \delta)$.

Series expansion for the symmetric case $\beta=0$

The starting point is the integral representation $$ F(x; \alpha, 0, \delta, \mu) = \frac{\delta e^{\delta \alpha}}{\sqrt{2\pi}}\int_0^{\infty} \Phi\left(\frac{x-\mu}{\sqrt{t}} \right)t^{-3/2} e^{-\frac{\alpha^2}{2}t} e^{-\frac{\delta^2}{2t}} \mathop{dt}. $$

We proceed expanding the term $e^{-\delta^2 / (2t)}$, $$ e^{-\frac{\delta^2}{2t}} = \sum_{k=0}^{\infty} \frac{(-1)^k}{2^k k!}\frac{\delta^{2k}}{t^k}. $$ Substituting this expansion into the integral and interchanging the order of summation and integration gives $$ F(x; \alpha, 0, \delta, \mu) = \frac{\delta e^{\delta \alpha}}{\sqrt{2\pi}} \sum_{k=0}^{\infty} \frac{(-1)^k \delta^{2k}}{2^k k!} I_k, $$ where $$ I_k = \int_0^{\infty} \Phi\left(\frac{x-\mu}{\sqrt{t}} \right)t^{-3/2-k} e^{-\frac{\alpha^2}{2}t}\mathop{dt}. $$

The integral $I_k$ converges when $x-\mu<0$. We therefore assume this condition and use the reflection identity otherwise, $$ F(x; \alpha, \beta, \delta, \mu) = 1 - F(-x; \alpha, -\beta, \delta, -\mu). $$

For $k\geq 1$, integration by parts gives $$ \begin{aligned} I_k &= -\frac{\alpha^2}{2k + 1}\int_0^{\infty}\Phi\left(\frac{x-\mu}{\sqrt{t}} \right)t^{-k-1/2} e^{-\frac{\alpha^2}{2}t}\mathop{dt} - \frac{x-\mu}{2k + 1}\frac{1}{\sqrt{2\pi}}\int_0^{\infty} e^{-\frac{(x-\mu)^2}{2t}} t^{-k-2} e^{-\frac{\alpha^2}{2}t}\mathop{dt}\\ &= -\frac{\alpha^2}{2k + 1} I_{k-1} - \frac{x-\mu}{2k + 1} \sqrt{\frac{2}{\pi}} \left(-\frac{\alpha}{x-\mu}\right)^{k+1} K_{k+1}(-\alpha (x-\mu)). \end{aligned} $$

Equivalently, $$ I_k = -\frac{\alpha}{2k+1}\left(\alpha I_{k-1} - \sqrt{\frac{2}{\pi}}\left(\frac{\alpha}{|x-\mu|}\right)^k K_{k+1}(\alpha|x-\mu|) \right). $$

The recurrence mainly involves computing the modified Bessel functions $K_0(x)$ and $K_1(x)$, followed by the Bessel recurrence to obtain the subsequent functions $K_{k+1}(x)$. The main difficulty is how to efficiently compute $I_0$ without resorting to numerical integration. We investigate that in the following section.

Computing $I_0$

Let $\xi=x-\mu<0$. Differentiating $I_k$ with respect to $\xi$ gives $$ \frac{\partial I_k}{\partial \xi} = \int_0^\infty \frac{1}{\sqrt{t}} \phi\left(\frac{\xi}{\sqrt{t}}\right) t^{-3/2-k}e^{-\alpha^2t/2}\mathop{dt}. $$

Since $\xi<0$, we write the result of evaluating the integral in terms of the positive argument $-\alpha\xi$: $$ \frac{\partial I_k}{\partial \xi} = \sqrt{\frac{2}{\pi}} \left(-\frac{\alpha}{\xi}\right)^{k+1} K_{k+1}(-\alpha\xi). $$

Because $I_k(\xi)\to0$ as $\xi\to-\infty$, integration with respect to $\xi$ yields $$ I_k = \sqrt{\frac{2}{\pi}} \int_{-\infty}^{\xi} \left(-\frac{\alpha}{t}\right)^{k+1} K_{k+1}(-\alpha t)\mathop{dt}. $$

Now apply the change of variables $y=-\alpha t$. Since $-\alpha t>0$ over the integration interval, this gives $$ I_k = \sqrt{\frac{2}{\pi}}\alpha^{2k+1} \int_{-\alpha\xi}^{\infty} \frac{K_{k+1}(y)}{y^{k+1}}\mathop{dy}. $$

For $k=0$, this reduces to $$ I_0 = \sqrt{\frac{2}{\pi}}\alpha \int_{-\alpha\xi}^{\infty} \frac{K_1(y)}{y}\mathop{dy}. $$

An antiderivative of $K_1(y)/y$ can be expressed in terms of modified Struve functions (Rosenheinrich, 2003): $$ \int \frac{K_1(y)}{y}\mathop{dy} = -yK_0(y)-K_1(y) -\frac{\pi y}{2} \left[ K_0(y)\mathbf{L}_1(y) + K_1(y)\mathbf{L}_0(y) \right], $$ where $\mathbf{L}_n(y)$ denotes a modified Struve function. Evaluating this antiderivative at the integration limits gives $$ I_0 = \sqrt{\frac{2}{\pi}}\alpha \left[ rK_0(r)+K_1(r) +\frac{\pi r}{2} \left( K_0(r)\mathbf{L}_1(r) + K_1(r)\mathbf{L}_0(r) \right) -\frac{\pi}{2} \right], $$ where $r=-\alpha\xi>0$.

Convergence analysis

To analyze the convergence of the series, we study the behavior of $I_k$ as $k\to\infty$. We start from the integral representation $$ I_k = \sqrt{\frac{2}{\pi}} \alpha^{2k+1} \int_{-\xi \alpha}^{\infty} \frac{K_{k+1}(y)}{y^{k+1}} \mathop{dy}. $$ For fixed $y>0$, the modified Bessel function has the following large-order asymptotic estimate $$ K_{k+1}(y) \sim \sqrt{\frac{\pi}{2(k+1)}} \left(\frac{e y}{2(k+1)} \right)^{-k-1}. $$

Substituting this approximation into the integral and taking absolute values gives $$ |I_k| \sim 2^{k+1}(k+1)^{k+1/2}e^{-k-1} \alpha^{2k+1} \int_{-\alpha\xi}^{\infty}y^{-2(k+1)}\mathop{dy} = \frac{2^{k+1}(k+1)^{k+1/2}e^{-k-1}} {(2k+1)|\xi|^{2k+1}}. $$

Consequently, the absolute values of the coefficients $$ c_k=\frac{\delta^{2k}}{2^k k!}I_k $$ satisfy $$ |c_k| \sim \left(\frac{\delta}{|\xi|}\right)^{2k} \frac{1}{2k+1} \left(\frac{k+1}{k}\right)^{k+1/2} \sqrt{\frac{2}{\pi}}\frac{1}{e|\xi|}. $$

Since $$ \lim_{k\to\infty} \left(\frac{k+1}{k}\right)^{k+1/2} =e, $$ we obtain $$ |c_k| \sim \frac{1}{2k+1} \sqrt{\frac{2}{\pi}}\frac{1}{|\xi|} \left(\frac{\delta}{|\xi|}\right)^{2k}. $$

Therefore, the series converges absolutely when $$ \delta<|\xi|=|x-\mu|. $$

Upper bound for the remainder

Let define the series $S = S_N + R_N$, where $$ S_N = \sum_{k=0}^{N-1} \frac{(-1)^k\delta^{2k}}{2^k k!}I_k $$ denote the truncated inner series, and let $$ R_N = \sum_{k=N}^{\infty} \frac{(-1)^k\delta^{2k}}{2^k k!}I_k $$ be its remainder. To discuss the size of the terms, define $$ a_k=\frac{\delta^{2k}}{2^k k!}I_k. $$ The sequence $I_k$ increases with $k$, whereas the prefactor $$ b_k=\frac{\delta^{2k}}{2^k k!} $$ is unimodal, since $$ \frac{b_{k+1}}{b_k}=\frac{\delta^2}{2(k+1)}. $$ Thus, the terms $a_k=b_kI_k$ need not decrease from the first term. However, the asymptotic estimate obtained above gives $$ a_k \sim \frac{1}{2k+1} \sqrt{\frac{2}{\pi}}\frac{1}{|x-\mu|} \left(\frac{\delta}{|x-\mu|}\right)^{2k}, $$ and therefore $$ \frac{a_{k+1}}{a_k} \longrightarrow \left(\frac{\delta}{|x-\mu|}\right)^2. $$ Therefore, in the convergent regime $\delta<|x-\mu|$, the term magnitudes are eventually decreasing. For truncation indices $N$ beyond this eventual-monotonicity threshold, the alternating-series remainder estimate gives $$ |R_N| \leq a_N = \frac{\delta^{2N}}{2^N N!}I_N. $$

For $x-\mu<0$, we use the standard normal-tail bound $$ \Phi(-z) \leq \frac{e^{-z^2/2}}{z\sqrt{2\pi}}, \qquad z>0. $$ This gives $$ \begin{aligned} I_N &= \int_0^{\infty} \Phi\left(\frac{x-\mu}{\sqrt{t}}\right) t^{-3/2-N}e^{-\alpha^2t/2},dt \\ &\leq \frac{1}{|x-\mu|} \int_0^{\infty} t^{-N-1} \exp\left( -\frac{(x-\mu)^2}{2t} -\frac{\alpha^2t}{2} \right),dt \\ &= \frac{2}{|x-\mu|} K_N\left(\alpha|x-\mu|\right) \left(\frac{|x-\mu|}{\alpha}\right)^N. \end{aligned} $$

Thus, $$ |R_N| \leq \frac{\delta^{2N}}{2^N N!} \frac{2}{|x-\mu|} K_N\left(\alpha|x-\mu|\right) \left(\frac{|x-\mu|}{\alpha}\right)^N. $$

A simpler, though less sharp, estimate follows by applying an upper bound to the modified Bessel function: $$ |R_N| \leq \left(\frac{\delta}{|x-\mu|}\right)^{2N} \frac{\alpha|x-\mu|}{2N}. $$

Numerical tests

The following examples illustrate the effectiveness of the proposed series expansion in the small-$\delta$ regime. We compare the number of terms required by the convergence-acceleration expansion described in Navas-Palencia (2026), Section 3.1.2 with those required by the proposed expansion. The reported errors are measured against a reference value computed with sufficiently high precision using mpmath.

Comparison with convergence acceleration for $|x-\mu|\to0$

The table shows that the proposed expansion requires fewer terms in the regime of study, the reduction being more noticeable for $\alpha \geq 1$. Given the simplicity of the series compared to the exponentially improved expansion, this should be preferred option for implementation. It is, however, interesting to remark that for very large $\alpha$, the recurrence for $I_k$ contains terms of magnitude proportional to $\alpha^2$. As a result, the computation may involve some cancellation, and the accuracy of the initial value $I_0$ becomes especially important.

The original Bessel-type expansion in Navas-Palencia (2026), equation (3.3.1) is absolutely convergent in the general case, but it can converge slowly when $\delta$ is small. The expansion presented here provides an alternative for this regime and can avoid the need for numerical integration. In practical implementations, even for small $\delta$, this expansion should be combined with asymptotic and uniform-asymptotic expansions.

The extension to the general case $\beta\neq0$ does not introduce substantial additional difficulties in the recurrence. However, the initial integral $I_0$ no longer appears to have a closed-form expression of the type obtained above.

Applications

The proposed expansion is useful whenever the scale parameter $\delta$ is small relative to the distance from the location parameter, that is, when $$ \delta<|x-\mu|. $$ This situation occurs naturally in short-maturity applications of the exponential NIG model, where the scale parameter is proportional to $\delta\tau$ and $\tau\ll1$. In this regime, the expansion can reduce the need for numerical integration and provide an efficient alternative for evaluating the NIG distribution function. See also the discussion of the exponential NIG model.

Small values of $\delta$ can also arise when NIG models are calibrated to financial time series, including energy-market data. In such applications, the expansion may improve the efficiency of repeated distribution-function evaluations within calibration, likelihood estimation, and pricing procedures. The use of the NIG distribution for energy-market modelling is discussed by Benth and Šaltytė-Benth (2004).

References

• Benth, F. E., and Šaltytė-Benth, J. (2004). The Normal Inverse Gaussian distribution and spot price modelling in energy markets. International Journal of Theoretical and Applied Finance, 7(2), 177–192.
• Navas-Palencia, G. (2026). On the computation of the cumulative distribution function of the normal inverse Gaussian distribution. Numerical Algorithms 102, 2349–2389.
• Rosenheinrich, W. (2003). Tables of some indefinite integrals of Bessel functions of integer order.