function liki = kalman(para,y);
sig = para(1);
bet = .99;
kap = para(2);
rho = para(3);
chi_pi = para(4);
chi_x = para(5);
rho_u = para(6);
rho_g = para(7);
sig_e = para(8);
sig_u = para(9);
sig_g = para(10);
gain = para(11);
xi_pi = (1-rho)*chi_pi;
xi_x = (1-rho)*chi_x;
T=length(y);
% The state is [pi,x,i,u,g,pi+1,x+1]
% The observable is [pi,x,i];
GAM0 = eye(7);
GAM0(1:2,:) = [1,-kap,0,-1,0,-bet,0;0,1,sig,0,-1,-sig,-1];
GAM0(6:7,6:7) = zeros(2,2);
GAM0(6:7,1:2) = eye(2);
GAM1 = zeros(7,7);
GAM1(3,:) = [xi_pi,xi_x,rho,0,0,0,0];
GAM1(4:7,4:7) = diag([rho_u,rho_g,1,1]);
PSI = zeros(7,3);
PSI(3:5,:) = eye(3);
PPI = zeros(7,2);
PPI(6:7,:) = eye(2);
Q_t = diag([sig_e^2,sig_u^2,sig_g^2]);
c = zeros(7,1);
[T_t,c_t,R_t,fmat,fwt,ywt,gev,eu,loose]=gensys(GAM0,GAM1,c,PSI,PPI,1+1e-10);
% function [G1,C,impact,fmat,fwt,ywt,gev,eu,loose]=gensys(g0,g1,c,psi,pi,div)
% System given as
%        g0*y(t)=g1*y(t-1)+c+psi*z(t)+pi*eta(t),
% with z an exogenous variable process and eta being endogenously determined
% one-step-ahead expectational errors.  Returned system is
%       y(t)=G1*y(t-1)+C+impact*z(t)+ywt*inv(I-fmat*inv(L))*fwt*z(t+1) .
P_t = inv(eye(49)-kron(T_t,T_t))*kron(R_t,R_t)*reshape(Q_t,9,1);
P_t = reshape(P_t,7,7);
phi = zeros(4,3,T);
R = zeros(4,4,T);
a_lt = zeros(3,1,T);
b_lt = zeros(3,3,T);
b_lt(:,:,1) = T_t(1:3,1:3);
phi(:,:,1) = [a_lt(:,:,1),b_lt(:,:,1)]';
R(1,1,1) = 1;
R(2:4,2:4,1) = P_t(1:3,1:3);
e_pi = zeros(T,1);
e_x = zeros(T,1);
e_z = zeros(3,1,T);
c_t = zeros(5,1,T);
for i = 2:T
    x = [1;y(:,i-1)];
    phi(:,:,i) = phi(:,:,i-1) + gain*inv(R(:,:,i-1))*x*(y(:,i)-phi(:,:,i-1)'*x)';
%     if rcond(phi(2:4,:,i))<1e-10
%         liki = -1e+6;
%         return
%     end
%     if rcond(phi(2:4,:,i))==NaN
%         liki = -1e+6;
%         return
%     end
    R(:,:,i) = R(:,:,i-1) + gain*(x*x'-R(:,:,i-1));
    a_lt(:,i) = phi(1,:,i)';
    b_lt(:,:,i) = phi(2:4,:,i)';
    if max(abs(eigs(b_lt(:,:,i))))>1
        phi(:,:,i) = phi(:,:,i-1);
        a_lt(:,i) = a_lt(:,i-1);
        b_lt(:,i) = b_lt(:,i-1);
    end
    e_z(:,i) = a_lt(:,i-1) + b_lt(:,:,i-1)*a_lt(:,i-1) + (b_lt(:,:,i-1))^2*y(:,i-1);
    e_pi(i) = e_z(1,i);
    e_x(i) = e_z(2,i);
    c_t(:,i) = [bet*e_pi(i);e_x(i)+sig*e_x(i);0;0;0];
end

if sum(eu) ~= 2
    liki = -1e+6;
    return
end
% Under a projection facility, models with adaptive learning generally
% imply that beliefs, if not outright converge to the REE, will form a
% stationary distribution around the REE. For this reason we will discard
% parameter draws which do not have a unique and stable REE
% Now we integrate learning into the model
% s_t = [pi,x,i,u,g]
% And we have a time-varying c_t that includes expectations, which are the
% heart of this model
Z_t = [eye(3),zeros(3,2)];
gam0 = eye(5);
gam0(1:2,:) = [1,-kap,0,-1,0;0,1,sig,0,-1];
gam1 = zeros(5,5);
gam1(3,:) = [xi_pi,xi_x,rho,0,0];
gam1(4:5,4:5) = diag([rho_u,rho_g]);
r_t = [zeros(2,3);eye(3)];
gam_inv = inv(gam0);
T_t = gam_inv*gam1;
R_t = gam_inv*r_t;
P_t = P_t(1:5,1:5);
a_t = zeros(5,1);
liki = 0;
% Remember that without using kalman-filtered states as in
% Slobodyan+Wouters, beliefs are pinned down by observables and the gain
% parameter. 
d_t = zeros(3,1);
H_t = zeros(3,3);
for i = 1:T
    a_t1 = T_t*a_t + c_t(:,i);
    P_t1 = T_t*P_t*T_t' + R_t*Q_t*R_t';
    y_t1 = Z_t*a_t1+d_t;
    v_t = y(:,i)-y_t1;
    F_t = Z_t*P_t1*Z_t' + H_t;
%     if rcond(F_t)<1e-10
%         liki = -1e+6;
%         return
%     end
    F_inv = inv(F_t);
    inv_F = F_inv;
    %update
    a_t = a_t1 + P_t1*Z_t'*inv_F*v_t;
    P_t = P_t1 - P_t1*Z_t'*inv_F*Z_t*P_t1;
    liki = liki -(3/2)*log(2*pi) - (1/2)*log(det(F_t)) - (1/2)*v_t'*F_inv*v_t;
end
if liki ~=real(liki)
    liki = -1e+6;
    return
end
