function llf = Likeli_Spatial_Panel(W,XC,Y,iter,ns,param)


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



% latent process speciified as: lambda_t = kappa*lambda_t-1 + theta*W*lambda_t-1 + rho*W*lambda_t + X_t*beta + tau + e_t    with   tau~N(0,sigma_tau) and e_t~N(0,sigma_e)
% param denotes Kx1 vector of paramters with  param(1)=kappa, param(2)=theta, param(3)=rho, param(4:K-2) = beta, param(K-1)=sigma_tau and param(K)=sigma_e
% W denotes NxN matrix of spatial weights -> should be declared as a sparse matrix!!
% XC cell array of regressors with T matrices of dimension N x K-5 
% Y denotes NxT matrix of a dependent variable (T is time dimension) 
% iter denotes number of EIS iterations 
% ns  denotes number of trajectories 



% create common random numbers (crn's)
rng(1234567) % control seed
N=length(W);
T=size(Y,2);
crnu=randn(ns,N*T)';


% obtain permutation vector ppp from precision matrix of eta=[lambda; tau]
ppp=permutation(N,T,W);


cholcheck=0; % later needed to check whether cholesky factor can be computed


% construct joint density of latent process with precision matrix "H" (see appendix 2 in the paper, here we develop everything in terms of sparse matrices!!)
par.kappa=param(1); par.theta=param(2); par.rho=param(3); par.beta=param(4:end-2); par.sigtau=param(end-1); par.sigmu=param(end);
d=1e-27; % => d is lower bound for sigma_tau and sigma_e
par.sigtau(par.sigtau<=d)=d; par.sigmu(par.sigmu<=d)=d;
Dinv=speye(N,N)-par.rho*W;
A=par.kappa*speye(N)+par.theta*W;
% stab_eig=max(abs(eig(full(Dinv\A)))); % process is stable if stab_eig<1
Ha=(1/par.sigmu^2)*(Dinv'*Dinv);


b=zeros(N*T-N,1);
% btau=spalloc(N*T-N,N,ceil(0.12*N^2*T));
% H=spalloc(N*T-N,N*T-N,ceil(0.12*(N*T)^2));
il=1; iu=N;
y=Y(:,1);
y(y==0)=0.5;
b(il:iu)=(1/par.sigmu^2)*Dinv'*(A*log(y) + XC{2}*par.beta);
btau(il:iu,1:N)=(1/par.sigmu^2)*Dinv';
atau=-(1/par.sigmu^2)*(A*log(y) + XC{2}*par.beta);
H(il:iu,il:iu)=Ha;
msta=A*log(y) + XC{2}*par.beta;
sm=msta'*(1/par.sigmu^2)*msta;

for t=3:T
    il=(t-2)*N+1; iu=(t-1)*N;
    b(il-N:iu-N)=-(A'/par.sigmu^2)*(XC{t}*par.beta) + b(il-N:iu-N);
    b(il:iu)=(1/par.sigmu^2)*Dinv'*(XC{t}*par.beta);
    btau(il-N:iu-N,1:N)=-(A'/par.sigmu^2) + btau(il-N:iu-N,1:N);
    btau(il:iu,1:N)=(1/par.sigmu^2)*Dinv';
    atau=atau-(1/par.sigmu^2)*XC{t}*par.beta;
    H(il-N:iu-N,il-N:iu-N)=H(il-N:iu-N,il-N:iu-N)+A'*(1/par.sigmu^2)*A;
    Hc=-Dinv'*(1/par.sigmu^2)*A;
    H(il:iu,il:iu)=Ha; H(il:iu,il-N:iu-N)=Hc; H(il-N:iu-N,il:iu)=Hc';
    sm=sm+par.beta'*XC{t}'*(1/par.sigmu^2)*XC{t}*par.beta;
end
GM=((T-1)/par.sigmu^2+(1/par.sigtau^2));
K=spdiags(GM*ones(N,1),0,N,N);
pppa=symamd(Ha);
Ha=Ha(pppa,pppa);
try
ldH = (T-1)*2*sum(log(diag(chol(Ha)))); % log determinant of Ha for t=2->T
catch
ldH = (T-1)*2*sum(log(diag(chol(spdiags(ones(N,1),0,N,N)))));
cholcheck=1;
end
s=sm+N*(T-1)*log(2*pi)-ldH;
y=reshape(Y,N*T,1);
y(1:N)=[];
m=H\b;


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


% P-matrix and mu-vector
Pmat = [(H + B) -btau; -btau' K];
qvec = [beta + b; atau];
mu = Pmat\qvec;


% reorder to reduce number of non-zero elements during EIS procedure
Pmat=Pmat(ppp,ppp);
mu=mu(ppp);


% simulating from the initial sampler
try
Lmat = chol(Pmat); % log determinant of Ha for t=2->T
catch
Lmat = chol(spdiags(ones(N*T,1),0,N*T,N*T));
cholcheck=1;
end
eta = repmat(mu,1,ns) + Lmat\crnu; % N x ns
% re-reorder and extract lambda
eta(ppp,:)=eta;
la=eta(1:N*(T-1),:);


% loop through the EIS iterations
 eisc = zeros(N*(T-1),3);
 for k=1:iter
     la = la'; % ns X N
     
     % EIS regressions
     for i=1:N*(T-1)
         yy = -gammaln(y(i)+1) + y(i)*la(:,i) - exp(la(:,i)); % ns x 1
         xx = [ones(ns,1) la(:,i) (la(:,i).*la(:,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*(T-1),N*(T-1));     
     beta = eisc(:,2);
     kappa = eisc(:,3);
     
     % P-matrix and q-vector
     Pmat = [(H + B) -btau; -btau' K];
     qvec = [beta + b; atau];
     mu = Pmat\qvec;
     
     
     % reorder to reduce number of non-zero elements during EIS procedure
     Pmat=Pmat(ppp,ppp);
     mu=mu(ppp);

     
     % simulating from the EIS sampler
     try
         Lmat = chol(Pmat); % log determinant of Ha for t=2->T
     catch
         Lmat = chol(spdiags(ones(N*T,1),0,N*T,N*T));
         cholcheck=1;
     end
     eta = repmat(mu,1,ns) + Lmat\crnu; % N x ns
     % re-reorder and extract lambda
     eta(ppp,:)=eta;
     la=eta(1:N*(T-1),:);
 end
 
 % computing the EIS estimate of the likelihood function
 ldPmat = 2*sum(log(diag(Lmat))); % log determinant of Pmat
 r = sum(kappa) + s + N*log(2*pi) - N*log(1/par.sigtau^2);
 
 
 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 - (N*(T-1)+N)*log(2*pi)   );
 fy = sum(-repmat(gammaln(y+1),1,ns) + repmat(y,1,ns).*la - exp(la)); 
 lng  =  0.5*( laB - 2*laC + repmat(sum(kappa),1,ns)) + fy; % 1 x ns
 scal =  mean(lng);
 lngs = lng - scal;
 tg1  = log(mean(exp(lngs)));
 llf  = -(scal + tg1 + lchi);

if isnan(llf)
    llf=1e50;
elseif isinf(llf)
    llf=1e50;
elseif isreal(llf)==0
    llf=1e50;
elseif llf<=0
    llf=1e50;
elseif cholcheck==1
    llf=1e50;
% elseif stab_eig>=1
%     llf=1e50;    
end

end



function ppp=permutation(N,T,W)

par.kappa=0.1; par.rho=0.1; par.theta=0.1; par.sigtau=0.2; par.sigmu=0.3;
Dinv=speye(N,N)-par.rho*W;
A=par.kappa*speye(N)+par.theta*W;

% btau=spalloc(N*T-N,N,ceil(0.12*N^2*T));
% H=spalloc(N*T-N,N*T-N,ceil(0.12*(N*T)^2));
il=1; iu=N;
btau(il:iu,1:N)=(1/par.sigmu^2)*Dinv';
Ha=(1/par.sigmu^2)*(Dinv'*Dinv);
H(il:iu,il:iu)=Ha;

for t=3:T
    il=(t-2)*N+1; iu=(t-1)*N;
    btau(il-N:iu-N,1:N)=-(A'/par.sigmu^2) + btau(il-N:iu-N,1:N);
    btau(il:iu,1:N)=(1/par.sigmu^2)*Dinv';
    H(il-N:iu-N,il-N:iu-N)=H(il-N:iu-N,il-N:iu-N)+A'*(1/par.sigmu^2)*A;
    Hc=-Dinv'*(1/par.sigmu^2)*A;
    H(il:iu,il:iu)=Ha; H(il:iu,il-N:iu-N)=Hc; H(il-N:iu-N,il:iu)=Hc';
end
GM=((T-1)/par.sigmu^2+(1/par.sigtau^2));
K=spdiags(GM*ones(N,1),0,N,N);
B=spdiags(ones(N*(T-1),1),0,N*(T-1),N*(T-1));
Pmat = [(H + B) -btau; -btau' K];
ppp=symamd(Pmat);

clear B DINV GM H HA HC K LA Pmat ans btau bu il iu Dinv Ha Hc

end