This paper deals with the estimation of normal inverse Gaussian distribution parameters from the Bayesian approach using Markov chain Monte Carlo methods. The normal inverse Gaussian distribution has the probability density function \[ f(y;\alpha,\beta,\mu,\delta)= {\alpha\over\pi}\exp(\delta\gamma+\beta(y-\mu))q(y)^{-1/2} K_{1}(\delta\alpha q(y)^{1/2}), \] where \(\alpha=\sqrt{\gamma^2+\beta^2}\), \(q(y)=1+[(y-\mu)/\delta]^2\), and \(K_{r}(x)\) denotes the modified Bessel function of the third kind of order \(r\). The authors use the data augmentation offered by the mixture derivation of the normal inverse Gaussian distribution which leads to Gibbs sampling scheme. Different prior distributions are examined, starting with a parsimonious conjugate prior in the augmented context. The basic normal inverse Gaussian model is extended to include covariates. Such a model is quite general and, among others, allows modelling various types of heteroscedastic and autoregressive processes. The authors apply the proposed Markov chain Monte Carlo scheme to a data set concerning the monthly nominal arithmetic returns of the FT-Actuaries All-Share Index for the UK.