function [Pit, Pic, Vt, Vc, Phi1, Phi2, V_ptr] = MCMC_inf(infl,PTR, parsinit, prior2, nobs, M)

% use PTR as an proxy to help identify inflation trend
% PTR = \pi* + e_ptr    e_ptr ~ N(0, V_ptr)

disp('#3) MCMC Inflation sampler...')
pars2 = parsinit ;

Pit = nan(M,nobs);
Pic = nan(M,nobs);
Vt = nan(M,1);
Vc = nan(M,1);
Phi1 = nan(M,1);
Phi2 = nan(M,1);
V_ptr = nan(M,1);


tic ;  % counting elapsed time 
for j = 1:M
    if mod(j,1000)==0
    fprintf('Iteration %d \n',j) ;
    toc ;  % duration so far
    tic ;
    end
    
    
    % 1) draw states -- pit and pic
    sim2 = KF2(infl,PTR, pars2, 1);
    pit = sim2(1,:)';
%     piy = sim2(2,2:nobs)';
%     pix = [sim2(3,1:nobs-1)' sim2(3,2:nobs)'];
    piy = sim2(2,2:nobs)';
%     pix = [sim2(2,2:nobs-1)' sim2(2,1:nobs-2)'];
    pix = [sim2(2,1:nobs-1)' sim2(3,1:nobs-1)'];
    
    % 2) phi1 phi2 ~ N
    phi_cov2 = (prior2.Omega_inv + pars2(2)\(pix'*pix)) \ eye(2) ; 
    [pchol, Sflag2] = chol(phi_cov2);
        if Sflag2 ~= 0 
           [UP, SP, ~] = svd(phi_cov2);
           pchol = (UP*sqrt(SP))';
        end
    phi_mean2 = pchol' * pchol * (prior2.Omega_inv * prior2.Phi + pars2(2)\(pix'*piy)) ;
    phi_draw2 = phi_mean2 + pchol' * randn(2,1) ;
    if abs(phi_draw2(2))<1 && sum(phi_draw2)<1 && phi_draw2(2)-phi_draw2(1)<1
        pars2(3:4) = phi_draw2 ;
    end
    
    % 3) vc ~ IG
    pic_hat = pix * pars2(3:4); 
    qerrors2 = piy -  pic_hat;
    pars2(2) = rinvgamma(1, (prior2.alpha_e2 + nobs-1)/2, (prior2.delta_e2 + sum(qerrors2.^2,'all'))/2 );
%     pars2(2) = rinvgamma(1, (prior2.alpha_e2 + nobs-2)/2, (prior2.delta_e2 + sum(qerrors2.^2,'all'))/2 );
    
    % 4) vt ~ IG
    qerrors1 = diff(pit);
    pars2(1) = rinvgamma(1, (prior2.alpha_e1 + nobs-1)/2, (prior2.delta_e1 + sum(qerrors1.^2,'all'))/2 );
    
    % 5) v_ptr ~ IG
    qerrors3 = PTR - pit;
    qerrors3 = qerrors3(~isnan(qerrors3));
    pars2(5) = rinvgamma(1, (prior2.alpha_e3 + size(qerrors3,1))/2, (prior2.delta_e3 + sum(qerrors3.^2,'all'))/2 );

    % save draws
    Pit(j,:) = pit;
    Pic(j,:) = sim2(2,:);
    Vt(j,1) = pars2(1);
    Vc(j,1) = pars2(2);
    V_ptr(j,1) = pars2(5);
    Phi1(j,1) = pars2(3);
    Phi2(j,1) = pars2(4);
     
end