function y = KF2(infl, PTR, pars, index)
% Kalman Filter to update beta_t|t and to generate log likelihood
% Two observed variables: inflation, PTR
% PTR is quarterly available, so has many nan

T = size(infl,1);
beta1 = zeros(T,3) ;    % beta t| t
beta2 = zeros(T,3) ;    % beta t| t-1
ptt1  = zeros(3*T,3) ;    % p   t | t
ptt2  = zeros(3*T,3) ;    % p   t | t-1
% parameterization
vt = pars(1) ;   % variance of the innovation to pit
vc = pars(2) ;   % variance of the innovation to pic
ph1 = pars(3);
ph2 = pars(4);
v_ptr = pars(5); % variance of the innovation to PTR

F = [1 0 0; 0 ph1 ph2; 0 1 0];   % transition matrix
F1 = F(2:3, 2:3);   % pick out AR part
H = [1 1 0; 1 0 0];        % measurement matrix
Q = zeros(3);       % variance-covariance of innovations in transit equation
Q(1,1) = vt;
Q(2,2) = vc;
Q1 = Q(2:3, 2:3);    % covariance of AR part
R = [0 0; 0 v_ptr]; 

% beta_ll = [mean(infl); 0; 0];   % starting values
% if ind ==1 
%     beta_ll = [3; 0; 0];   % starting values
% elseif ind ==2
%     beta_ll = [6; 0; 0];   
% elseif ind==3
%     beta_ll = [15; 0; 0];   
% end
beta00 = mean(infl(1:12));
beta_ll = [beta00; infl(1)-beta00; 0]; 

% variacne matrix of initial state variables
p_ll = zeros(3,3);
p_ll1 = (eye(4)-kron(F1,F1)) \ reshape(Q1,4,1);
p_ll1 = reshape(p_ll1', 2,2);
p_ll(1,1) = 100;
p_ll(2:3,2:3) = p_ll1;
j_iter = 1;

while j_iter <= T
    if isnan(PTR(j_iter))
        beta_tl = F*beta_ll;
        p_tl = F *p_ll *F' + Q;
        vt = infl(j_iter,1) - H(1,:)*beta_tl;         % prediction error for inflation
        ft = H(1,:)*p_tl*H(1,:)' + 0;                 % variacne of forecast error
        beta_tt= beta_tl + (p_tl * H(1,:)' * vt)/ft;
        p_tt= p_tl - (p_tl * H(1,:)' * H(1,:) * p_tl)/ft;
    else
        beta_tl = F*beta_ll;
        p_tl = F *p_ll *F' + Q;
        vt = [infl(j_iter,1); PTR(j_iter)] - H*beta_tl;    % prediction error
        ft = H*p_tl*H' + R;                                % variacne of forecast error
        beta_tt= beta_tl + p_tl * H' /(ft)* vt;
        p_tt= p_tl - p_tl * H' /(ft)* H * p_tl;
    end

    beta1(j_iter,:) = [beta_tt(1,1) beta_tt(2,1) beta_tt(3,1)] ;
    beta2(j_iter,:) = [beta_tl(1,1) beta_tl(2,1) beta_tl(3,1)] ;
    ptt1(3*(j_iter-1)+1:3*j_iter,:) = p_tt ;
    ptt2(3*(j_iter-1)+1:3*j_iter,:) = p_tl ;


    beta_ll = beta_tt;
    p_ll = p_tt;

    j_iter = j_iter +1;
    
end
% llf = -sum(lik_mat);

if index==0
    y = -llf ;
else
% generating sequence of state variables (beta1 beta2,...,betaT)
% based on Kim and Nelson eq (8.9), (8.16'), (8.17') when Q is singular
%     stt = beta1(2:end,:)';  % 3x(T-1)
%     pttb = ptt1(4:end,:);   % 3x(3(T-1))
    stt = beta1';  % 3xT
    pttb = ptt1;   % 3x3T
    sTT = beta1(end,:)';
    pTT = ptt1(end-2:end,:);
    Fs = F(1:2,:);
    Qs = Q(1:2,1:2);
    
    [pchol, Sflag2] = chol(pTT);
    if Sflag2 ~= 0
        [UP, SP, ~] = svd(pTT);
        pchol = (UP*sqrt(SP))';
    end
    sT = sTT + pchol'*randn(3,1);
    s = zeros(3,T);
    s(:,T) = sT(1:3);
    
    for jj = 1:T-1
        sTm1 = stt(:,T-jj); %sTm1=st|t
        pTm1 = pttb(3*(T-jj-1)+1:3*(T-jj),:); %pTm1=pt|t
        J = pTm1*Fs'/(Fs*pTm1*Fs'+Qs);
        sTm1T = sTm1+J*(sT(1:2)-Fs*sTm1);     %sTm1T=st|t,st+1
        pTm1T = pTm1-J*Fs*pTm1;               %pTm1T=pt|t,st+1
        [pchol, Sflag2] = chol(pTm1T);
        if Sflag2 ~= 0 
           [UP, SP, ~] = svd(pTm1T);
           pchol = (UP*sqrt(SP))';
        end
        sT = sTm1T+pchol'*randn(3,1);
        s(:,T-jj) = sT;
    end
    y = s;
end