
function [beta_mcmc] = drawB(y_t,x_t,SIGMA_t,M,K,p,t,lambda,beta_0_prmean,beta_0_prvar)

beta_pred = zeros(K,t);
beta_update = zeros(K,t);
C_ttm1 = zeros(K,K,t);
C_tt = zeros(K,K,t);
e_t = zeros(M,t);

% 1\ Kalman filter  
for irep = 1:t    
    % Predict
    if irep==1
        beta_pred(:,irep) = beta_0_prmean;
        C_ttm1(:,:,irep) = beta_0_prvar;
    else
        beta_pred(:,irep) = beta_update(:,irep-1);
        C_ttm1(:,:,irep) = (1./lambda)*C_tt(:,:,irep-1); %C_tt(:,:,irep-1) + OMEGA;          
    end
    
    e_t(:,irep) = y_t(irep,:)' - x_t((irep-1)*M+1:irep*M,:)*beta_pred(:,irep);
   
    %update beta[t]
    Rx = C_ttm1(:,:,irep)*x_t((irep-1)*M+1:irep*M,:)';
    KV = SIGMA_t(:,:,irep) + x_t((irep-1)*M+1:irep*M,:)*Rx;
    KG = Rx/KV;       
    beta_update(:,irep) = beta_pred(:,irep) + KG*e_t(:,irep); %#ok<*MINV>
    C_tt(:,:,irep) = C_ttm1(:,:,irep) - KG*(x_t((irep-1)*M+1:irep*M,:)*C_ttm1(:,:,irep));            
end

% 2/ Rauch-Tung-Striebel fixed-interval smoother for the time-varying parameters  
beta_new = 0*beta_update; beta_new(:,t) = beta_update(:,t);  beta_mcmc = 0*beta_update;
C_tt_new = 0*C_tt; C_tt_new(:,:,t) = C_tt(:,:,t);
beta_mcmc(:,t) = mvnrnd(beta_new(:,t),C_tt_new(:,:,t),1);
for irep = t-1:-1:1
    % BETA_t
    U_t = lambda;%squeeze(C_tt(:,:,irep)/C_ttm1(:,:,irep+1));      
%     beta_new(:,irep) = beta_update(:,irep) + U_t*(beta_mcmc(:,irep+1) - beta_pred(:,irep+1));       
%     C_tt_new(:,:,irep) = C_tt(:,:,irep) - U_t*C_tt(:,:,irep);
    beta_new(:,irep) = beta_update(:,irep) + U_t*(beta_new(:,irep+1) - beta_pred(:,irep+1));       
    C_tt_new(:,:,irep) = C_tt(:,:,irep) + U_t*(C_tt(:,:,irep+1) - C_ttm1(:,:,irep+1))*U_t';     
    beta_mcmc(:,irep) = mvnrnd(beta_new(:,irep),C_tt_new(:,:,irep),1);
            
    % Check stationarity of draw
    bbtemp = beta_mcmc(M+1:K,irep);  % get the draw of B(t) at time i=1,...,T  (exclude intercept)
    By = zeros(M,M*p);   
    for ii = 1:p              
        By(:,(ii-1)*M+1:ii*M) = reshape(bbtemp((ii-1)*(M^2)+1:ii*(M^2),1),M,M)';
    end
    B = [By ; eye(M*(p-1)) zeros(M*(p-1),M)];
    % Accept-Reject sampler for draw
    while max(abs(eig(B)))>0.999
        beta_mcmc(:,irep) = mvnrnd(beta_new(:,irep),C_tt_new(:,:,irep),1);       
        bbtemp = beta_mcmc(M+1:K,irep);  % get the draw of B(t) at time i=1,...,T  (exclude intercept)
        By = zeros(M,M*p);   
        for ii = 1:p              
            By(:,(ii-1)*M+1:ii*M) = reshape(bbtemp((ii-1)*(M^2)+1:ii*(M^2),1),M,M)';                
        end
        B = [By ; eye(M*(p-1)) zeros(M*(p-1),M)];
    end    
end