simulations - Use NIG distribution to model stock path

I would like to use Monte Carlo simulation to price some options. First I use standard approach where stock price is discribed by the following process: $$S_T = S_0\exp \left[(r - 0.5\sigma^2)T + \sigma \sqrt{T}\varepsilon\right],$$ where $S_0$ is an initial stock price at time $t = 0$, $S_T$ is a stock price at time $T$, $T$ is a time step, $r$ is the risk free interest rate, $\sigma$ is a standart deviation if stock returns and $\varepsilon$ is an independent normally distributed variable $\varepsilon \sim \phi(0, 1)$ There is an assumption under this model that log returns of a stock price have a normal distribution, such that $\log{\frac{S_T}{S_0}} \sim \phi(\mu, \sigma^2)$, where (from Ito's Lemma) $\mu = (r - 0.5\sigma^2)T$ and $var = \sigma^2T . $

Now I assume that returns of a stock price follow Normal Inverse Gaussian distribution (NIG). My question is am I able just to plug into my formula for $S_T$ NIG proccess with estimated parameters (using MLE) instead of $\varepsilon$ and than compare option prices? Proccess will look as following: $$S_T = S_0\exp \left[(r - 0.5\sigma^2)T + \sigma T\varepsilon^*\right],$$ where $\varepsilon^* \sim NIG(\hat{\mu}, \hat{\alpha}, \hat{\delta}, \hat{\beta}).$

I guess that properties of NIG don't allow me to do such manipulations but I am not sure where to find the solution to this problem. Any hints are welcome!

---

Many thanks @Quantuple for help:

Answer

Assuming you've used this definition for the NIG distribution and that you've managed to come up with estimates $(\hat{\alpha}, \hat{\beta}, \hat{\mu}, \hat{\delta} )$ of the individual NIG parameters, your question boils down to:

"How to simulate paths from the global log-return process $R_t = \ln(S_t/S_0)$ for all $t \in [0,T]$, assuming i.i.d. $NIG(\hat{\alpha}, \hat{\beta}, \hat{\mu}, \hat{\delta} )-$distributed periodic log-returns (in your case daily)?"

First of all, no, you cannot use the equation you mention.

Yet, because the NIG distribution is a special case of normal variance-mean mixture (see this document, page 14), if one lets \begin{align} \sigma^2 &\sim IG\left( \frac{\delta}{\gamma}, \delta^2 \right),\ \ \text{with } \gamma = \sqrt{ \alpha^2 - \beta^2 } \\ \varepsilon &\sim \mathcal{N}(0,1) \end{align} then the random variable $X$ defined as $$ X = \mu + \beta \sigma^2 + \sigma \varepsilon $$ follows a $NIG(\alpha,\beta,\mu,\delta)$ distribution.

Now, let $r_{\delta t, i}$ denote the $\delta t$-period log-return observed for a given $t_i \in [0,T]$ $$ r_{\delta t, i} := \ln\left( \frac{S_{t_i}}{S_{t_i-\delta t}} \right) $$

You can then proceed as follows to generate realisations of the global return process $(R_t)_{t\geq 0}$.

1. 
   

   Under the assumptions of i.i.d. NIG-distributed $\{r_{\delta t, i}\}$, estimate the NIG parameters $(\hat{\alpha}, \hat{\beta}, \hat{\mu}, \hat{\delta} )$ from historical data using your favourite method (Maximum Likelihood Estimation, Moment Matching etc.).

   
   


2. 
   

   To simulate the periodic log-returns $\{r_{\delta t,i}\}_{i=1,...,N}$ (in practice $N = T/\delta t$ where $T$ figures the horizon of your MC simulation and $\delta t$ the length of the period over which the individual periodic log-returns prevail) use the key result mentioned above by (i) generating $\sigma_i^2 \sim IG(\hat{\delta}/\hat{\gamma}, \hat{\delta}^2)\ \ \text{i.i.d.}$ (see here for instance), (ii) generating $\varepsilon_i \sim \mathcal{N}(0,1) \ \ \text{i.i.d.}$ (iii) computing $$ r_{\delta t,i} = \hat{\mu} + \hat{\beta} \sigma_i^2 + \sigma_i \varepsilon_i $$

   
   


3. 
   

   Once all $\{r_{\delta t, i}\}_{i=1,...,N}$ have been simulated, build the global return process $(R_t)_{t\geq 0}$ by aggregating a certain number $n$ of periodic returns. Indeed, for any fixed $t \in [0,T] $, the global log-return $R_t := \ln(S_t/S_0)$ computes as: \begin{align} R_t &= \ln\left(\frac{S_t}{S_0}\right) \\ &= \ln\left(\frac{S_t}{S_{t-\delta t}} \frac{S_{t-\delta t}}{S_{t- 2\delta t}} \dots \frac{S_{\delta t}}{S_{0}}\right) \\ &= \sum_{ t_i \in [\delta t, t] } \ln \left( \frac{S_{t_i}}{S_{t_i-\delta t}} \right) \\ &= \sum_{i \leq n} r_{\delta t,i} \end{align}

   
   
   

This method was first proposed by Rydberg in:

T. H. Rydberg. The normal inverse Gaussian Lévy process: simulation and approximation. Comm. Statist. Stochastic Models, 13(4):887–910, 1997.