function llf= Likeli_EIS_Poisson_Code(W,X,y,iter,ns,par)


% W denotes NxN matrix of spatial weight => W should be declared as a sparse matrix!! 
% y denotes Nx1 vector of a dependent count variable 
% iter denotes number of EIS regressions 
% ns  denotes number of trajectories 
% par denotes Kx1 vector of paramters with par(1) being the spatial 
... paramter rho and par(end) being the standard deviation of the latent process!!

% the output argument llf denotes the log-likelhood estimate


warning('off', 'MATLAB:rankDeficientMatrix')
warning('off', 'MATLAB:nearlySingularMatrix')


% determine number of observations
N=length(W);

% starting values
param=par;
% applies sine restriction to spatial paramter rho which ensures |rho|<=c
c=1-0.00001;
param(1)=rhorestrict(par(1),c); % => search value for param(1) which satisfies sin(param(1))*c=rho

% create common random numbers (crn's)
rng(1234567) % control seed
crnu=randn(ns,N)';


% expected value "m" and inverse covariance matrix "Smat" for SAL model
Destinv=speye(N,N)-sin(param(1))*c*W;
Smat=(1/param(end)^2)*(Destinv'*Destinv);
m=Destinv\(X*param(2:end-1));
% reorder to reduce number of non-zero elements during EIS procedure
ppp=symamd(Smat);
Smat=Smat(ppp,ppp);
m=m(ppp);
y=y(ppp);

        

% initial sampler based upon  TSE approximations
alpha=exp(m);
beta=exp(m).*(m-1);
B=spdiags(alpha,0,N,N);

% P-matrix and mu-vector
Pmat = Smat + B;
mu = Pmat\(y + beta - m.*alpha);


% simulating from the initial sampler
Lmat = chol(Pmat);
u = repmat(mu,1,ns) + Lmat\crnu;


 % loop through the EIS iterations
 eisc = zeros(N,3);
 for k=1:iter
     u = u';
     omm = repmat(m',ns,1) + u; % ns x N
     
     % EIS regressions
     for i=1:N         
         yy = - exp(omm(:,i)); % ns x 1
         xx = [ones(ns,1) omm(:,i) (omm(:,i).*omm(:,i))]; % ns x 3
         beis = xx\yy;
         beis(isnan(beis))=inf; % needed since NaN doesn't allow chol(Pmat)
         alpha = -2*beis(3); 
         alpha(alpha<0)=0; % assures positive definiteness of Pmat
         beta  =  beis(2);
         kappa = -2*beis(1);
         
         % store EIS parameters needed for likelihood evaluation
         eisc(i,1) = alpha;
         eisc(i,2) = beta;
         eisc(i,3) = kappa;
     end
     
     alpha = eisc(:,1);
     B=spdiags(alpha,0,N,N);     
     beta = eisc(:,2);
     ka = eisc(:,3);
     
     % P-matrix and q-vector
     Pmat = Smat + B;
     mu = Pmat\(y + beta - m.*alpha);
     
     % simulating from the EIS sampler
     Lmat = chol(Pmat);
     u = repmat(mu,1,ns) + Lmat\crnu;
 end
 
 % computing the EIS estimate of the likelihood function
 la = repmat(m,1,ns) + u;
 ldPmat = 2*sum(log(diag(Lmat))); % log determinant of Pmat
 Lmat = chol(Smat);
 ldSmat = 2*sum(log(diag(Lmat))); % log determinant of Smat
 r = sum(alpha.*(m.^2) -2*beta.*m + ka - 2*y.*m + 2*gammaln(y+1));
 
 laB = zeros(1,ns);
 laC = zeros(1,ns);
 for s=1:ns
     laB(s) = la(:,s)'*B*la(:,s);
     laC(s) = la(:,s)'*beta;
 end
 
 lchi =  -0.5*(r - mu'*Pmat*mu + ldPmat - ldSmat   );
 lng  =  0.5*( laB - 2*laC + repmat(sum(ka),1,ns)) - sum(exp(la)); % 1 x ns
 scal =  mean(lng);
 lngs = lng - scal;
 tg1  = log(mean(exp(lngs)));
 llf  = -(scal + tg1 + lchi); % negative log-likelihood estimate
end


function x=rhorestrict(rho,c)

 x= fsolve(@sinerestrict,0.3,optimset('fsolve'),rho,c);
 
end
 

function yr=sinerestrict(x,rho,c)
    
    yr=rho-sin(x)*c;

end