function [X,z] = nigrnd(N, alpha, beta, delta, mu,seed)
%the parameters should be vector of the same length
%N = number of simulated values
%X = column vector length(N) with ith entry a simulated
%    value from the NIG dist with parameters: alpha,
%    beta, delta and mu 
% z is the sensitivity dX/dbeta
% written by Joe Dicesare
rand('state',seed);
X = zeros(1,N);
one1 = ones(1,N);
chi = delta^2;
phi = alpha^2 - beta^2;
nu = sqrt(chi/phi);
v = chi2rnd(1,1, N);
z1 = nu*one1 + (nu^2/2/chi)*v  - nu/2/chi*sqrt(4*nu*chi*v+ nu^2*v.^2);
U = rand(1,N);
p = nu*(nu*one1+z1).^(-1);
H = (U <= p);
z = z1.*H + (nu^2*z1.^(-1)).*(one1 - H);
randn('state',seed);
y = randn(1,N);
X = mu*one1 +  beta*z + sqrt(z).*y;  % z  has the inverse gaussian distribution