function [SIGMA_t] = drawSIGMA(y_t,x_t,beta_mcmc,M,t,delta,SIGMA_0,S_0,n_0)

S_t = zeros(M,M,t);
n_t = zeros(1,t);
SIGMA_t = zeros(M,M,t);

e_t = zeros(M,t);

% 2.1. Forward recursions  
for irep = 1:t
    e_t(:,irep) = y_t(irep,:)' - x_t((irep-1)*M+1:irep*M,:)*beta_mcmc(:,irep);  
    e2_t = e_t(:,irep)*e_t(:,irep)';
    if irep == 1
        n_t(irep) = delta*n_0 + 1;  % Degrees of freedom
        a_t = 1/n_t(irep);
        S_t(:,:,irep) = (1-a_t)*S_0 + a_t*( (S_0^(0.5))*(SIGMA_0^(-.5))*e2_t*(SIGMA_0^(-.5))*(S_0^(0.5)) );                                    
    else
        n_t(irep) = delta*n_t(irep-1) + 1;  % Degrees of freedom
        a_t = 1/n_t(irep);
        S1 = sqrtm(S_t(:,:,irep-1));
        S2 = (SIGMA_t(:,:,irep-1)^(-.5)); % sqrtm(inv(SIGMA_t(:,:,irep-1)))
        S_t(:,:,irep) = (1-a_t)*S_t(:,:,irep-1) + a_t*( S1*S2*e2_t*S2*S1 );   
    end
    SIGMA_t(:,:,irep) = symmetric(S_t(:,:,irep))/(n_t(irep)-M-1);%iwishrnd(symmetric(S_t(:,:,irep)),n_t(irep));
end

% 2.2. Backward recursions   
S_t_new = 0*S_t;    S_t_new(:,:,t) = S_t(:,:,t);
n_t_new = 0*n_t;    n_t_new(:,t) = n_t(:,t);  
SIGMA_t(:,:,t) = iwishrnd(symmetric(S_t_new(:,:,t)),n_t_new(t));
for irep = t-1:-1:1       
    n_t_new(:,irep) = (1-delta)*n_t(:,irep) + delta*n_t_new(:,irep+1);
    S_t_new(:,:,irep) = (1-delta)*S_t(:,:,irep) + delta*S_t_new(:,:,irep+1);
    SIGMA_t(:,:,irep) = iwishrnd(symmetric(S_t_new(:,:,irep)),n_t_new(irep));
end