function [kinfQ, lamQ, gamma, Xbar, sigtau2, sige2, Phi, Omega, Sigma_tilde, tau, X, Alpha, Pit, Pic, Vt, Vc, phi1, phi2, Vptr, RMSE] = MCMCsampler(Y, infl, PTR, parsinit, prior, mats, dt, rlb, nobs, M)
% MCMC sampler to draw parameters and state variables

disp('#3) MCMC sampler...')
pars = parsinit ;
N = 3 ;
nu = 5   ;    % for independence proposals, multivariate t draws
iter = 1000;
tol = 0.01;
scale_thetaQ = [100, 0.1, 1, 1]' ;
scale_ap = [1, 1, 0.1, 0.1]' ;

kinfQ= nan(M,1) ;
lamQ = nan(M,N) ;
gamma = nan(M,N) ;
Xbar = nan(M,N) ;
sigtau2 = nan(M,1) ;
sige2 = nan(M,1) ;
Phi = nan(M,N^2) ;
Omega = nan(M, N*(N+1)/2) ;
Sigma_tilde = nan(M, N*(N+1)/2) ;
tau = nan(M, nobs) ;
X = nan(M*N, size(Y,1)) ;
Alpha = zeros(M,5) ;  % Metropolis-Hastings acceptance probabilities
ltri = tril(ones(3)) ;

Pit = nan(M,nobs);
Pic = nan(M,nobs);
Vt = nan(M,1);
Vc = nan(M,1);
phi1 = nan(M,1);
phi2 = nan(M,1);
Vptr = nan(M,1);

RMSE = nan(M,length(mats));

thetaQ_H0 = 0.01*eye(4);
sigma_H0 = 0.01*eye(6);
sigtau2_H0 = 0.1^2;
ap_H0 = 0.01*eye(4);

rng(123)
for j = 1:M
    % 1) draw states -- tau(i-star) and X
    sim = EKF(pars, Y, rlb, mats, pars.aJ, pars.bJ, pars.sigmasJ, 3);
    pars.tau = sim(1,:)' ;
    pars.X = sim(2:end,:)' ;  % tx3
    pars.llkQ = llk_Q(Y, pars.X, pars.aJ, pars.bJ, pars.sigmasJ, pars.sige2, rlb) ;
    pars.llkP = llk_P(pars.tau, pars.X, pars.Xbar, pars.Phi, pars.gamma, pars.omega, pars.sigtau2, nobs) ;
    % draw states -- pit and pic
    sim2 = KF2(infl, PTR, [pars.vt;pars.vc;pars.phi1;pars.phi2; pars.vptr], 1);
    pit = sim2(1,:)';
    piy = sim2(2,2:nobs)';
    pix = [sim2(2,1:nobs-1)' sim2(3,1:nobs-1)'];
    RMSE(j,:) = rmse_y(Y, pars.X, pars.aJ, pars.bJ, pars.sigmasJ, pars.sige2, rlb) ;
    
    % 2) thetaQ --- kinfQ + lamQ -- IMH
    lprior_current = log(normpdf(pars.kinfQ, prior.kmean, prior.kvar));
    thetaQ_current = [pars.kinfQ; log(-pars.dlamQ./(1+pars.dlamQ))].* scale_thetaQ; 
    [~,respars,~,resh,~,~,~] = csminwel(@(parameters)obj_thetaQ(parameters,Y,pars,prior,mats,dt,rlb,scale_thetaQ),thetaQ_current,thetaQ_H0,[],tol,iter);   
    thetaQ_mean = respars ;
    thetaQ_mvtscale = makePD(resh);
    corr_thetaQ = corrcov(thetaQ_mvtscale);
    thetaQ_prop = thetaQ_mean + 0.001*mvtrnd(corr_thetaQ, nu)' ;
    prop = prop_thetaQ(thetaQ_prop, Y, pars, prior, mats, dt, rlb, scale_thetaQ) ;
    if isstruct(prop)
        lr_prop = log(mvtpdf(thetaQ_current - thetaQ_mean, corr_thetaQ, nu)) - ...
            log(mvtpdf(thetaQ_prop - thetaQ_mean, corr_thetaQ, nu)) ;
        lr_prior = prop.lprior - lprior_current;
        prop.lprior = nan ;
        Alpha(j,1) = min(exp(prop.llkQ - pars.llkQ + lr_prop + lr_prior), 1) ;
        if rand(1)<Alpha(j,1)
            pars.kinfQ = prop.kinfQ;
            pars.lamQ = prop.lamQ;
            pars.dlamQ = prop.dlamQ;
            pars.llkQ = prop.llkQ ;
            pars.aJ = prop.an ;
            pars.bJ = prop.bn ;
            pars.sigmasJ = prop.SigmaQ ;
            thetaQ_H0 = thetaQ_mvtscale;
        end
    end

    % 3) Sigma -- Independence MH
    sigma_current = chol(pars.omega,'lower') ;
    sigma_current = sigma_current(ltri==1) ;
    lprior_current = log(diwish(pars.Sigma_tilde, prior.kappa, prior.Psi));
    [~,respars,~,resh,~,~,~] = csminwel(@(parameters)obj_sigma(parameters,Y,pars,prior,mats,dt,rlb,nobs), sigma_current, sigma_H0,[],tol,iter);
    sigma_mean = respars ;
    sigma_mvtscale = makePD(resh);
    corr_sigma = corrcov(sigma_mvtscale);
    sigma_prop = sigma_mean + 0.005*mvtrnd(corrcov(sigma_mvtscale), nu)';
    prop = prop_sigma(sigma_prop, Y, pars, prior, mats, dt, rlb, nobs) ;
    if isstruct(prop)
        lr_prop = log(mvtpdf(sigma_current - sigma_mean, corr_sigma, nu)) - ...
            log(mvtpdf(sigma_prop - sigma_mean, corr_sigma, nu))  ;
        lr_prior = prop.lprior - lprior_current ;
        prop.lprior = nan ;
        Alpha(j,2) = min(exp(prop.llkQ - pars.llkQ + prop.llkP - pars.llkP + ...
            lr_prop + lr_prior), 1) ;
        if rand(1)<Alpha(j,2)
            pars.omega = prop.omega ;
            pars.Sigma_tilde = prop.Sigma_tilde ;
            pars.llkQ = prop.llkQ ;
            pars.llkP = prop.llkP ;
            pars.aJ = prop.an ;
            pars.sigmasJ = prop.SigmaQ ;
            sigma_H0 = sigma_mvtscale;
        end
    end
    
    % 4) ap (gamma/Xbar) -- Independence MH  
    lprior_current = log(diwish(pars.Sigma_tilde, prior.kappa, prior.Psi)) + ...
        log(mvnpdf([pars.a; pars.p], prior.ap_mean, prior.ap_cov)) ;
    ap_current = [pars.a ; pars.p];
    [~,respars,~,resh,~,~,~] = csminwel(@(parameters)obj_ap(parameters,pars,prior,nobs, scale_ap), ap_current, ap_H0,[],tol,iter);
    ap_mean = respars ;
    ap_cov = makePD(resh) ;
    ap_prop = ap_mean +  10*chol(ap_cov,'lower') * randn(4,1) ;  % 4x1
    prop = prop_ap(ap_prop, pars,  prior, nobs, scale_ap) ;
    if isstruct(prop)
        lr_prop = logrationorm(ap_current, ap_prop, ap_mean, ap_cov) ;
        lr_prior = prop.lprior - lprior_current ;
        prop.lprior = NaN ;
        Alpha(j,3) = min(exp(prop.llkP - pars.llkP + lr_prop + lr_prior),1) ;
    end
    if rand(1)<Alpha(j,3)
        pars.a = prop.a ;
        pars.p = prop.p ;
        pars.gamma = prop.gamma ;
        pars.Xbar = prop.Xbar ;
        pars.Sigma_tilde = prop.Sigma_tilde ;
        pars.llkP = prop.llkP ;
        ap_H0 = ap_cov;
    end
    
    % 5) sigtau2 -- Independence MH
    lprior_current = log(dinvgamma(pars.sigtau2, prior.alpha_eta/2, prior.delta_eta/2)) + ...
        log(diwish(pars.Sigma_tilde, prior.kappa, prior.Psi)) ;
    [~,respars,~,resh,~,~,~] = csminwel(@(parameters)obj_sigtau2(parameters,pars,prior,nobs), log(pars.sigtau2), sigtau2_H0,[],tol,iter);
    lsigtau2_mean = respars ;
    lsigtau2_var = resh;
    lsigtau2_prop = lsigtau2_mean + 4*sqrt(lsigtau2_var)*trnd(nu,1) ;
    lr_prop = log(tpdf((log(pars.sigtau2) - lsigtau2_mean)/sqrt(lsigtau2_var),nu)) - ...
        log(tpdf((lsigtau2_prop - lsigtau2_mean)/sqrt(lsigtau2_var),nu)) ;
    sigtau2_prop = exp(lsigtau2_prop) ;
    prop = prop_sigtau2(sigtau2_prop, pars, prior, nobs)  ;
    if isstruct(prop)
        lr_prior = prop.lprior - lprior_current ;
        prop.lprior = NaN ;
        Alpha(j,4) = min(exp(prop.llkP - pars.llkP + lr_prior + lr_prop), 1) ;
        if rand(1)<Alpha(j,4)
            pars.sigtau2 = prop.sigtau2 ; 
            pars.Sigma_tilde = prop.Sigma_tilde ;
            pars.llkP = prop.llkP ;
            sigtau2_H0 = resh;
        end
    end
    
    % 6) Phi -- Gibbs sampling
    Xtilde = pars.X -  repmat(pars.Xbar', nobs,1) - pars.tau .* pars.gamma' ;  % Tx3
    Ymat = Xtilde(4:nobs,:)' ; % 3x(T-2)
    Zmat = Xtilde(3:nobs-1,:)' ;
    Sigma_inv = inv(pars.Sigma_tilde) ;
    phi_cov = (prior.Omega_inv + kron((Zmat*Zmat'), Sigma_inv)) \ eye(N^2) ; 
    phi_mean = phi_cov * (prior.Omega_inv * reshape(prior.Phi,[],1) + ...
        kron(Zmat, Sigma_inv) * reshape(Ymat,[],1)) ;
    phi_draw = phi_mean + chol(phi_cov) * randn(N^2,1) ;
    Phi_draw = reshape(phi_draw, N, N) ;
    ev = max(abs(eig(Phi_draw))) ;
    Alpha(j,5) = ev ;
    if ev<1
        pars.Phi = Phi_draw ;
        pars.llkP = llk_P(pars.tau, pars.X, pars.Xbar, pars.Phi, pars.gamma, pars.omega, pars.sigtau2, nobs) ;
    end

    
    % 7) sige2 - measurement error variance
    %    Gibbs-sampling
    Yhat =  Gsr(pars.X, pars.aJ, pars.bJ, pars.sigmasJ, rlb) ;  % TxK
    Qerrors = Y - Yhat ;
    Qerrors = Qerrors(3:end);
    pars.sige2 = rinvgamma(1, (prior.alpha_e + nobs*length(mats))/2, (prior.delta_e + sum(Qerrors.^2,'all'))/2 ) ;
    pars.llkQ = -0.5*sum(Qerrors.^2,'all')/pars.sige2 ;
    
    % 8) phi1 phi2 ~ N
    % Gibbs-sampling
    phi_cov2 = (prior.Omega2_inv + pars.vc\(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 * (prior.Omega2_inv * prior.Phi2 + pars.vc\(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
        pars.phi1 = phi_draw2(1) ;
        pars.phi2 = phi_draw2(2) ;
    end
    
    % 9) vc ~ IG
    % Gibbs-sampling
    pic_hat = pix * [pars.phi1; pars.phi2]; 
    qerrors2 = piy -  pic_hat;
    pars.vc = rinvgamma(1, (prior.alpha_e2 + nobs-1)/2, (prior.delta_e2 + sum(qerrors2.^2,'all'))/2 );
    
    % 10) vt ~ IG
    % Gibbs-sampling
    qerrors1 = diff(pit);
    pars.vt = rinvgamma(1, (prior.alpha_e1 + nobs-1)/2, (prior.delta_e1 + sum(qerrors1.^2,'all'))/2 );

    % 5) v_ptr ~ IG
    qerrors3 = PTR - pit;
    qerrors3 = qerrors3(~isnan(qerrors3));
    pars.vptr = rinvgamma(1, (prior.alpha_e3 + size(qerrors3,1))/2, (prior.delta_e3 + sum(qerrors3.^2,'all'))/2 );


    
    % save draws
    kinfQ(j) = pars.kinfQ ;
    lamQ(j,:) = pars.lamQ ;
    gamma(j,:) = pars.gamma ;
    Xbar(j,:) = pars.Xbar ;
    sigtau2(j) = pars.sigtau2 ;
    sige2(j) = pars.sige2 ;
    Phi(j,:) = reshape(pars.Phi, 1,[]) ;
    Omega(j,:) = pars.omega(ltri==1)  ;
    Sigma_tilde(j,:) = pars.Sigma_tilde(ltri==1)  ;
    tau(j,:) = pars.tau ;
    X((j-1)*N+1:j*N,:) = pars.X' ; 
    Pit(j,:) = pit;
    Pic(j,:) = sim2(2,:);
    Vt(j,1) = pars.vt;
    Vc(j,1) = pars.vc;
    phi1(j,1) = pars.phi1;
    phi2(j,1) = pars.phi2;
    Vptr(j,1) = pars.vptr;
end

end

























