The exponential NIG model for option pricing

In this post, we review the works in (Aguilar, J. P., 2020) and (Aguilar, J. P., & Kirkby, J. L., 2022) on the derivation of closed-form pricing formulas for a variety of path-independent options in the exponential normal inverse Gaussian (NIG) model. In these works, the authors obtain series expansions using the Mellin transform of the payoff and the NIG probability density function. In the following, we show that all path-independent options considered can be expressed in terms of the cumulative distribution function (CDF) of the NIG and the generalized hyperbolic (GH) distribution. These results can be easily implemented using statistical libraries that support the NIG and GH distributions (for example, SciPy stats.norminvgauss and stats.genhyperbolic, respectively), continuously leveraging any improvements in the algorithms for computing the corresponding CDFs.

For numerical tests, we reproduce the results in (Aguilar, J. P., 2020) using the C++ implementation of the algorithm developed in (Navas-Palencia, G., 2025) for the computation of the CDF of the NIG distribution. Furthermore, we devise a simple algorithm using two of the Bessel-type series derived in (Navas-Palencia, G., 2025) to obtain, at least for the range of parameters generally occurring in financial applications, accurate results up to single-machine precision (1e-8), an error acceptable in most cases.

The normal inverse Gaussian distribution

The NIG distribution is defined as a normal variance-mean mixture distribution, where the mixing distribution for the variance follows an inverse Gaussian (IG) distribution. Consider \(Z\) is a positive random variable following a distribution, \(Z \sim \mathcal{IG}(\delta^2, \gamma^2)\), \(\mu\) and \(\beta\) are constants. If \(X \sim \mathcal{NIG}(\alpha, \beta, \delta, \mu)\), then $$ X \,|\, Z \sim \mathcal{N}(\mu + \beta Z, Z), $$ where \(\gamma = \sqrt{\alpha^2 - \beta^2}\). The NIG distribution parameters admit the following interpretation: \(\alpha\) is the tail heaviness, \(\beta\) is the asymmetry or skewness, \(\mu\) is the location and \(\delta\) the scaling. If \(\beta = 0\), the distribution is symmetric and centred around \(\mu\). The domain of the parameters is $$ 0 \le |\beta| < \alpha, \quad \mu \in \mathbb{R}, \quad \delta > 0. $$ The NIG distribution has probability density function $$ f(x; \alpha, \beta, \delta, \mu) = \frac{\alpha \delta}{\pi} \frac{K_1\left(\alpha\sqrt{\delta^2 + (x-\mu)^2}\right)}{\sqrt{\delta^2 + (x-\mu)^2}} e^{\delta \sqrt{\alpha^2 - \beta^2} + \beta(x-\mu)}, $$ where \(K_1(x)\) is the modified Bessel function of the second kind. The cumulative distribution function has the common integral representation from its density function, $$ F(x; \alpha, \beta, \delta, \mu) = \frac{\alpha \delta e^{\delta \gamma}}{\pi} \int_{-\infty}^{x} \frac{K_1\left(\alpha\sqrt{\delta^2 + (t-\mu)^2}\right)}{\sqrt{\delta^2 + (t-\mu)^2}} e^{\beta(t-\mu)} \mathop{dt}. $$ On the other hand, using the variance-mean mixture distribution properties, we have $$ F(x; \alpha, \beta, \delta, \mu) = \frac{\delta e^{\delta\gamma}}{\sqrt{2\pi}}\int_{0}^{\infty} \Phi\left(\frac{x - (\mu +\beta t)}{\sqrt{t}}\right) t^{-3/2} e^{-\frac{\delta^2}{2t} - \frac{\gamma^2 t}{2}} \mathop{dt}, $$ where \(\Phi(x)\) is the standard normal cumulative distribution function. From the latter integral, it follows from the reflection formula of \(\Phi(x) = 1 - \Phi(-x)\) that \(F(x; \alpha, \beta, \mu, \delta) = 1- F(-x; \alpha, -\beta, -\mu, \delta)\).

The exponential NIG model

In the subsequent sections, we provide a summary of (Aguilar, J. P., 2020) closely following its expository style, and we refer to the paper for further details. Consider \(T > 0\) and \(S: t \in [0, T] \rightarrow S_t\) be the market price of some financial asset. The instantaneous variation of \(S_t\) is given by $$ \frac{\mathrm{d} S_t}{S_t} = (r- q) \mathrm{d} t + \mathrm{d} X_t. $$ where \(r \ge 0\) is the risk-free interest rate and \(q \ge 0\) is the dividend yield, and where \(\left\{X_t\right\}_{t \in [0, T]}\) is the NIG process. The solution to the stochastic differential equation is the exponential process $$ S_T = S_t e^{(r-q+\omega)\tau + X_{\tau}}, \quad \omega= -\mu + \delta\left(\sqrt{\alpha^2 - (\beta + 1)^2} - \sqrt{\alpha^2 - \beta^2} \right), $$ where \(\tau = T-t\) is the time to maturity, \(\omega\) is the martingale adjustment, and \(X_{\tau} \sim NIG(\alpha,\beta,\delta\tau, \mu\tau)\). Note that if \(Y_{\tau} = X_{\tau} - \mu \tau \sim NIG(\alpha,\beta,\delta\tau, 0)\), then the value of \(\mu\) is irrelevant for pricing, therefore, we set \(\mu = 0\).

Pricing contigent claims

Given a path-independent payoff function \(\mathcal{P}\) and some strike parameters \(K_1,\ldots, K_n > 0\), then the value at time \(t\) of a contingent claim delivering a payoff \(\mathcal{P}\) at maturity is equal to the following risk-neutral expectation: $$ \mathcal{C} = \mathbb{E}^{\mathbb{Q}}\left[e^{-r \tau} \mathcal{P}(S_T, K_1,\ldots,K_n) \,|\, S_t\right]. $$ The conditional expectation has the integral representation $$ \mathcal{C} = e^{-r\tau}\int_{-\infty}^{\infty} \mathcal{P}(S_t e^{(r-q+\omega)\tau + x}, K_1,\ldots,K_n) f(x; \alpha,\beta, \delta \tau, \mu \tau) \mathop{dx}. $$ where \(f(x; \alpha,\beta, \delta \tau, \mu \tau)\) is the probability density function of the NIG process.

Digital option (cash-or-nothing)

A cash-or-nothing call digital option is a type of binary option that pays a fixed amount of cash if the underlying asset \(S_T\) is above a specified strike price \(K\) at expiration, and nothing otherwise. The payoff is given by $$ \mathcal{P}_{c/n}(S_T, K) = \mathbb{1}_{\{S_T > K\}} $$ The value of a cash-or-nothing call is expressed in terms of the NIG CDF as follows $$ C_{c/n} = e^{-r\tau} F(k; \alpha, -\beta, \delta\tau, 0). $$ To derive the above expression, the first step in defining the log forward moneyness \(k\), calculated as

$$ S_T > K \Longleftrightarrow S_t e^{(r-q-\omega)\tau + X_{\tau}} > K \Longleftrightarrow X_{\tau} > -k, \quad k = \log{\frac{S_t}{K}} + (r-q + \omega)\tau. $$ then $$ \mathbb{E}[\mathbb{1}_{\{S_T > K\}}] = 1 - \mathbb{P}[X_t \le -k] = 1 - F(-k; \alpha, \beta, \delta\tau, 0) = F(k; \alpha, -\beta, \delta\tau, 0), $$ where the last step uses the reflection formula of the NIG CDF.

Digital option (asset-or-nothing)

An asset-or-nothing call digital option is a binary option that pays the value of the underlying asset \(S_T\) if it is above a specified strike price \(K\) at expiration, and nothing otherwise. The payoff can be written as

$$ \mathcal{P}_{a/n}(S_T, K) = S_T \mathbb{1}_{\{S_T > K\}} $$

The value of an asset-or-nothing call can be expressed in terms of the NIG CDF after replacing the expression below in its first integral representation presented above $$ \mathcal{P}_{a/n}\left(S_t e^{(r-q+\omega)\tau + x}, K\right) = K e^{k + x} \mathbb{1}_{\{x > -k\}}, $$ obtaining $$ C_{a/n} = S_t e^{-q\tau} F(k; \alpha, -\beta-1, \delta\tau, 0). $$

Note that the representation in terms of the NIG CDF imposes the constraint \(|\beta| + 1 < \alpha\). Nevertheless, Table 1 (Aguilar, J. P., 2020) illustrates that \(\alpha\) tends to be significantly larger than \(|\beta|\) for the implied NIG parameters obtained by calibrating various option markets.

European option

The European call option payoff is given by $$ \mathcal{P}_{eur}(S_T, K) = \left(S_T - K\right)^+, $$ and can be written as $$ \mathcal{P}_{eur}\left(S_t e^{(r-q+\omega)\tau + x}, K\right) = K \left(e^{k + x} - 1\right) \mathbb{1}_{\{x > -k\}}. $$ Therefore, the option price can be expressed in terms of the cash-or-nothing and the asset-or-nothing as follows $$ C_{eur} = C_{a/n} - K C_{c/n} = S_t e^{-q\tau} F(k; \alpha, -\beta-1, \delta\tau, 0) - K e^{-r\tau} F(k; \alpha, -\beta, \delta\tau, 0). $$

Log options

The procedure to get similar expressions for the log options is slightly more involved. The payoff of a log call is $$ \mathcal{P}_{log\;call}(S_T, K) = \left(\log S_T - \log K\right)^+, $$ and following the previous notation, we have $$ \mathcal{P}_{log\;call}\left(S_te^{(r - q + \omega)\tau + x}, K\right) = (k + x)^+. $$

Asymmetric model

First, we derive the pricing formulas for the asymmetric case \((\beta \neq 0)\). The conditional expectation integral can be decomposed as follows $$ C_{log} = e^{-r\tau}\int_{-\infty}^{\infty} (k+x)^+ f(x; \alpha, \beta, \delta\tau, 0) \mathop{dx} = k F(k; \alpha, -\beta, \delta\tau, 0) + \int_{-k}^{\infty} x f(x; \alpha, \beta, \delta\tau, 0) \mathop{dx}. $$ Now observe that the rightmost integral is expressible in terms of the NIG CDF and its derivative w.r.t. \(\beta\): $$ \int_{-k}^{\infty} x f(x; \alpha, \beta, \delta\tau, 0) \mathop{dx} = \frac{\partial F(k; \alpha, -\beta, \delta\tau, 0)}{\partial \beta} + \frac{\beta \delta\tau}{\gamma} F(k; \alpha, -\beta, \delta\tau, 0). $$ To calculate the derivative w.r.t. \(\beta\), we differentiate under the integral sign the variance-mean integral representation $$ \frac{\partial F(k; \alpha, -\beta, \delta\tau, 0)}{\partial \beta} = \frac{\delta\tau}{\sqrt{2\pi}} \int_0^{\infty} \frac{\partial}{\partial \beta}\left\{ \Phi\left(\frac{x + \beta t}{\sqrt{t}}\right) e^{-\frac{\gamma^2}{2}t + \delta\tau\gamma}\right\} t^{-3/2} e^{-\frac{(\delta\tau)^2}{2t}} \mathop{dt} = \frac{\delta\tau e^{\delta\tau\gamma}}{\sqrt{2\pi}}\left(I_1 + I_2\right), $$ where the integrals \(I_1\) and \(I_2\) are defined as $$ I_1 = \frac{\delta\tau e^{\delta\tau\gamma}}{\sqrt{2\pi}} \int_0^{\infty} \phi\left(\frac{k + \beta t}{\sqrt{t}}\right) t^{-1} e^{-\frac{\gamma^2}{2}t -\frac{(\delta\tau)^2}{2t}} \mathop{dt} $$ $$ I_2 = \frac{\delta\tau e^{\delta\tau\gamma}}{\sqrt{2\pi}} \int_0^{\infty} \Phi\left(\frac{k + \beta t}{\sqrt{t}}\right) \beta\left(t - \frac{\delta\tau}{\gamma}\right)t^{-3/2} e^{-\frac{\gamma^2}{2}t -\frac{(\delta\tau)^2}{2t}} \mathop{dt} $$ The integral \(I_1\) can be represented in terms of the modified Bessel function \(K_0(x)\) and \(I_2\) is a sum of a NIG CDF and GH CDF, the latter being the cumulative distribution function of the generalized hyperbolic distribution denoted as \(F_{GH}(x;\lambda,\alpha,\beta,\delta,\mu)\) with \(\lambda = 1/2\). After rearranging terms and simplifying, we get $$ C_{log} = e^{-r\tau} \left[k F(k;\alpha, -\beta, \delta \tau, 0) + \frac{\delta \tau e^{\gamma \delta\tau -\beta k}}{\pi} K_0\left(\alpha\sqrt{k^2 + (\delta\tau)^2}\right) + \frac{\beta \delta\tau}{\gamma}F_{GH}\left(k; \frac{1}{2}, \alpha, -\beta, \delta\tau, 0\right)\right]. $$

Symmetric model

Setting \(\beta = 0\), we remove the GH CDF, resulting in a much simpler expression $$ C_{log} = e^{-r\tau} \left[k F(k;\alpha, 0, \delta \tau, 0) + \frac{\delta \tau e^{\alpha \delta\tau}}{\pi} K_0\left(\alpha\sqrt{k^2 + (\delta\tau)^2}\right) \right]. $$

Power options

The power option payoff is a power function with exponent \(a > 0\) of the underlying asset's price \(S_T\) at expiration.

Cash-or-nothing

The payoff of the digial cash-or-nothing is $$ \mathcal{P}_{pow,\,c/n}(S_T, K) = \mathbb{1}_{\{S_T^a > K\}}. $$ As previously described, the corresponding log forward moneyness \(k_a\) is simply $$ k_a = \log\frac{S_t}{K^{1/a}} + (r-q + \omega)\tau $$ Proceeding as before, we obtain a cash-or-nothing call in terms of the NIG CDF. $$ C_{pow,\,c/n} = e^{-r\tau} F(k_a; \alpha, -\beta, \delta\tau, 0). $$

Asset-or-nothing

In a similar manner, we get the asset-or-nothing call $$ \mathcal{P}_{pow,\,a/n}(S_T, K) = S_T^a \mathbb{1}_{\{S_T^a > K\}} = Ke^{a(k_a + x)} \mathbb{1}_{\{x \ge -k_a\}} $$ $$ C_{pow,\,a/n} = K e^{-r\tau + a k_a - \delta\tau\left(\sqrt{\alpha^2 - (\beta + a)^2} - \sqrt{\alpha^2 - \beta^2}\right)} F(k_a; \alpha, -\beta-a, \delta \tau, 0), $$

and the representation in terms of the NIG CDF imposes the constraint \(|\beta| + a < \alpha\).

Implementation

The algorithms described in (Navas-Palencia, G., 2025) combine series expansions, uniform and asymptotic expansions and numerical integration, resulting in a robust implementation capable of attaining precision close to 1e-15. However, if the log forward moneyness \(k\) is small, two Bessel-type series can be used effectively to cover a wide range of the parameters' domain. For the special case \(\beta=0\), one can use the absolutely convergent series

$$ F(x; \alpha, 0, \delta, 0) = \frac{1}{2} + \frac{x\alpha \delta e^{\delta \alpha}} {\pi \sqrt{x^2 + \delta^2}}\sum_{j=0}^{\infty} \frac{1}{(2j +1)!!} \left(\frac{x^2 \alpha}{\sqrt{x^2 + \delta^2}}\right)^{j}K_{j+1}(\alpha \sqrt{x^2 + \delta^2}). $$

We note that a similar series expansion with alternating terms, although not absolutely convergent, is presented in (Aguilar, J. P., 2020) and (Navas-Palencia, G., 2025). For the general case, we use the series expansion

$$ F(x;\alpha, \beta, \delta, 0) = \frac{1}{2} + \frac{x\alpha\delta e^{\delta \gamma + x\beta}}{\pi \sqrt{x^2 + \delta^2}}\sum_{j=0}^{\infty} \frac{1}{(2j+1)!!}\left(\frac{x^2\alpha}{\sqrt{x^2 + \delta^2}}\right)^j A_j, $$ where $$ A_j = \sum_{m=0}^{2j+1} (-1)^m \binom{2j+1}{m} \left(\frac{\sqrt{x^2 + \delta^2} \beta}{\alpha x}\right)^m K_{j + 1 - m}(\alpha \sqrt{x^2 + \delta^2}). $$

Note that the terms of both series can be computed recursively by applying properties of the modified Bessel function. In addition, the convergence of the series improves for large values of \(\alpha\) and/or \(\delta\). A rigorous derivation of the remainder after series truncation and convergence analysis is discussed in (Navas-Palencia, G., 2025). Furthermore, other series expansions for small \(x\) and/or \(\beta\) in terms of Hermite polynomials are also derived, which might require fewer terms for some sets of parameters.

Numerical tests

To conclude, we replicate some of the numerical tests presented in (Aguilar, J. P., 2020). For these tests, we compare the reference implementation in (Navas-Palencia, G., 2025) with a simpler implementation combining the two Bessel-type expansions presented above. As demonstrated, fewer than 10 terms are sufficient to achieve a relative error below 1e-4, and fewer than 20 terms suffice for an error below 1e-8.

Symmetric: Table 4 \([\beta = 0]\) (Aguilar, J. P., 2020). Parameters: \(K=4000\), \(r=0.01\), \(q=0\), \(\tau=2\), \(\alpha=8.9932\) and \(\delta=1.1528\).

Asymmetric: Table 4 \([\beta = -4.5176]\) (Aguilar, J. P., 2020). Parameters: \(K=4000\), \(r=0.01\), \(q=0\), \(\tau=2\), \(\alpha=8.9932\) and \(\delta=1.1528\).

Asymmetric: Table 3 \([\beta = -4.5176]\) (Aguilar, J. P., 2020). Parameters: \(K=4000\), \(r=0.01\), \(q=0\), \(\tau=1\), \(\alpha=8.9932\) and \(\delta=1.1528\).

Finally, we reproduce the experiment shown in Figure 2 (Aguilar, J. P., 2020) for various maturities \(\tau \in \{1/12, 1/6, 1, 2\}\). The rest of parameters are \(S_t = 3000\), \(\alpha = 40\), \(\beta=0\), \(\delta = 25\), \(q=0\) and \(r=0.01\). The Python code can be found here.

The experiments are performed on an AMD Ryzen 7 5825U running at 2.0 GHz. The following table shows the time for computing 3000 European call options (\(K\in[2000, 5000]\)), comparing the Bessel-type expansion with the reference implementation. As discussed earlier, as \(\delta \tau \to \infty\), the Bessel-type expansions behave as asymptotic expansions, accelerating the convergence. For those cases, the CPU time for the call option is below 1 microsecond (600-900 nanoseconds). For smaller \(\delta\tau\), a slightly smaller working tolerance should be used to guarantee the desired user's tolerance.

The results illustrate that the reference implementation prioritizes accuracy, and for small \(\delta\tau\), numerical integration is chosen more often, being approximately 100 times slower. In this respect, the next release of the reference implementation will expose the desired user tolerance as an input and include additional uniform asymptotic expansions to reduce the regions where numerical integration is employed.

References

• Aguilar, J. P. (2021). Explicit option valuation in the exponential NIG model. Quantitative Finance, 21(8), 1281–1299.
• Aguilar, J. P., & Kirkby, J. L. (2022). Closed-form option pricing for exponential Lévy models: a residue approach. Quantitative Finance, 23(2), 251–278.
• Navas-Palencia, G. (2025). On the computation of the cumulative distribution function of the normal inverse Gaussian distribution. Numerical Algorithms.