I didn't look at your simulation results because IMHO you already have various problems on a strictly theoretical basis.
(1)
Writing $$ S_T = S_0 \exp \left[(r - 0.5\sigma^2)dt + \sigma dt\varepsilon\right] $$ makes little sense (mix of $T$ and $dt$ and wrong scaling of variance), I guess you meant: $$ S_{t+\delta t} = S_{t} \exp \left[(r - 0.5\sigma^2)\delta t + \sigma \sqrt{\delta t} \varepsilon\right] $$ which is - for instance - what is postulated in the Black-Scholes model (along with the independence of the $\{\epsilon\}$)
(2)
In that case, it is not $S_{t+\delta t}/S_t$ which is normally distributed but rather $\delta t$-periodic log-returns $r_{\delta t}$ $$ r_{\delta t} := \ln\left( \frac{S_{t+\delta t}}{S_t} \right) \sim \mathcal{N}(m,s^2) $$ \begin{align} m &= (r-\frac{\sigma^2}{2})\delta t \\ s^2 &= \sigma^2 \delta t \end{align}
(3)
If you estimate a sample mean and standard deviation ($\hat{m}$ and $\hat{s}$) from observed $\delta t$-period log-returns' time series (i.e. when working under the physical measure $\mathbb{P}$), you should better use the formulation $$ S_{t+\delta t} = S_{t} \exp \left[\hat{m} + \hat{s} \varepsilon\right] $$ with $\varepsilon \sim \mathcal{N}(0,1)$.
Your BS-like formulation indeed only makes real sense when working under a risk-neutral measure $\mathbb{Q}$.
(4)
You cannot write $$ S_{t+\delta t} = S_{t} \exp \left[\hat{m} + \hat{s} \varepsilon^* \right] $$ with $\varepsilon^*$ a random variable following a NIG distribution fitted to your data and reasonably expect that the $\delta t$-period log-returns follow NIG with mean $\hat{m}$ and variance $\hat{s}^2$ (because $\epsilon^*$ presumably does not exhibit zero mean and unit variance).
Still, assuming you have estimated your NIG parameters correctly (MLE or moment matching) you could write $$ r_{\delta t} \sim NIG(\hat{\alpha}, \hat{\beta}, \hat{\mu}, \hat{\delta} ) $$
(5)
Sampling realisations of a NIG process is described in equation (1) (page 6). Using different notations, you could write that $$ X = \mu + \beta \sigma^2 + \sigma \varepsilon $$ with \begin{align} \sigma^2 &\sim IG(\delta,\gamma) \\ \varepsilon &\sim \mathcal{N}(0,1) \end{align} and $$ \gamma = \sqrt{ \alpha^2 - \beta^2 } $$
follows a $NIG(\alpha,\beta,\mu,\delta)$ distribution.