function [L]=likeliehl(dy,Theta)

% Copyright: Pau Rabanal.  Feel free to copy, modify and use at your own risk.
% However, you are not allowed to sell this software or otherwise impinge
% on its free distribution.

% the following are the imput of the function

[T,Ny]=size(dy);

x0=[0;0;0;0;0;0;0;0;0]; % it has to be a way of calculating the steady-state of the state-variables
			 % given some Theta
Nz=length(x0);

D = zeros(Ny);

sigma  = Theta(1);%            %intertemporal elast. of substitution
thetap = Theta(2);%        %probability of not adjusting price
thetaw = Theta(3);%         %probability of not adjusting wages
gam    = Theta(4);%           %Coefficient on labor u. function
phi    = Theta(5);%       %elasticity of labor demand w.r.t. relative real wage
alpha  = Theta(6);%   		%rule of thumb parameter

% Coefficients on interest rate rule

rhor   = Theta(7);%  %Interest Rate Smoothing
gy     = Theta(8);%  %Coefficents in Taylor rule: output
gpi    = Theta(9);%  %								    price inflation
rhoa   = Theta(10);% % autocorrelation of technology shocks
rhog   = Theta(11);% % autocorrelation of demand shocks

% endogenous error. Note that this is specified in st.des. 

sigma_a= Theta(12);%  %tech. shock 
sigma_m= Theta(13);%  %mont. shock
sigma_p= Theta(14);%  %price shock 
sigma_g= Theta(15);%  %demand shock

Sigma = [sigma_a, 0, 0, 0,
         0, sigma_m, 0, 0,
         0, 0, sigma_p, 0,
         0, 0, 0, sigma_g];
% solving the model

[PP,QQ,RR,SS,NN,PROBLEM]=modelehl(sigma,thetap,thetaw,gam,phi,alpha,rhor,gy,gpi,rhoa,rhog);

if PROBLEM==1
        L=-Inf;
        return;
end
% writting the state-space form matrices
C=[QQ*Sigma; Sigma];
Ao=[PP QQ*NN; zeros(size(NN,1),size(PP,2)) NN];
G=zeros(Ny,Nz); %G selects dp,r,y from the state vector
G(1,1)=1;
G(2,2)=1;
G(3,3)=1;
G(4,5)=1;
%useful transformation
Gbar = G*Ao-D*G;
RR=zeros(4,4);
Rbar=RR+G*C*C'*G';
V = C*C';
% likelihood function
% intial values for the system
exog=size(Ao,1);
Sigma     =inv(eye((exog^2))-kron(Ao,Ao))*reshape(V,(exog^2),1);
Sigma     =reshape(Sigma,exog,exog);
xt        = x0';
innov     = dy(1,:)-xt*Gbar';
L          = 0;
for i=2:T;   
   Omega        = Rbar+Gbar*Sigma*Gbar';
   if det(Omega)==0;
       Omega=eye(4);
       L=-Inf;
       return;
   end
   Omegai       = inv(Omega);
   K            = (Ao*Sigma*Gbar'+C*C'*G')*Omegai;
   Sigma        = Ao*Sigma*Ao'+C*C'-K*Omega*K';
   L            = L + log(det(Omega))+innov*Omegai*innov';
   xt           = xt*Ao'+innov*K';
   innov        = dy(i,:)-xt*Gbar';
end
L              = L + log(det(Omega))+innov*Omegai*innov'+4*T*log(2*pi);
L=-(0.5*L);