function [Ht,h2,capAt,Sig2,Sblockdraw] = drawSIGMAmcmc(y,Z,Btdraw,h2,M,numa,p,t,Sig2,Sblockdraw,A_0_prmean,A_0_prvar,S_prmean,S_prvar,sigma_prmean,sigma_prvar,alpha2,beta2)


    % Substract from the data y(t), the mean Z x B(t)
    yhat = zeros(M,t);
    for i = 1:t
        yhat(:,i) = y(:,i) - Z((i-1)*M+1:i*M,:)*Btdraw(:,i);
    end
    
    % This part is more tricky, check Primiceri
    % Zc is a [M x M(M-1)/2] matrix defined in (A.2) page 845, Primiceri
    Zc = - yhat';
    sigma2temp = exp(h2);
    
    Atdraw = [];
    ind = 1;
    sizeS = 1:M;
    for ii = 2:M
        % Draw each block of A(t)
        [Atblockdraw,log_lik2a] = carter_kohn(yhat(ii,:),Zc(:,1:ii-1),sigma2temp(:,ii),...
            Sblockdraw{ii-1},sizeS(ii-1),1,t,A_0_prmean(((ii-1)+(ii-3)*(ii-2)/2):ind,:),A_0_prvar(((ii-1)+(ii-3)*(ii-2)/2):ind,((ii-1)+(ii-3)*(ii-2)/2):ind));
        Atdraw = [Atdraw ; Atblockdraw]; %#ok<AGROW> % Atdraw is the final matrix of draws of A(t)
        ind = ind + ii;
    end
    
    %=====| Draw S, the covariance of A(t) (from iWishart)
    % Take the SSE in the state equation of A(t)
    Attemp = Atdraw(:,2:t)' - Atdraw(:,1:t-1)';
    sse_2 = zeros(numa,numa);
    for i = 1:t-1
        sse_2 = sse_2 + Attemp(i,:)'*Attemp(i,:);
    end
    % ...and subsequently draw S, the covariance matrix of A(t) 
    ijc = 1;
    for jj=2:M
        Sinv = inv(sse_2(((jj-1)+(jj-3)*(jj-2)/2):ijc,((jj-1)+(jj-3)*(jj-2)/2):ijc) + S_prmean{jj-1});
        Sinvblockdraw = wish(Sinv,t+S_prvar(jj-1));
        Sblockdraw{jj-1} = inv(Sinvblockdraw); % this is a draw from S
        ijc = ijc + jj;
    end
    
    %------------------------------------------------------------------------------------------
    %   STEP III: Draw diagonal VAR covariance matrix log-SIGMA(t)
    %------------------------------------------------------------------------------------------
    % First create capAt, the lower-triangular matrix A(t) with ones on the
    % main diagonal. This is because the vector Atdraw has a draw of only the
    % non-zero and non-one elements of A(t) only.
    capAt = zeros(M*t,M);
    for i = 1:t
        capatemp = eye(M);
        aatemp = Atdraw(:,i);
        ic=1;
        for j = 2:M
            capatemp(j,1:j-1) = aatemp(ic:ic+j-2,1)';
            ic = ic + j - 1;
        end
        capAt((i-1)*M+1:i*M,:) = capatemp;
    end
    
    y2 = [];
    for i = 1:t
        ytemps = capAt((i-1)*M+1:i*M,:)*yhat(:,i);
        y2 = [y2  ((ytemps + 1e-10).^2)]; %#ok<AGROW>
    end
    
    Ytilde = log(y2)';
    for j=1:M
        [h2(:,j) Sg2(:,j)] = SVRW(Ytilde(:,j),h2(:,j),Sig2(j),sigma_prmean(j),sigma_prvar(j,j));
    end
    sigt = exp(.5*h2);

    e2 = h2(2:end,:) - h2(1:end-1,:);     
    S_1 = alpha2 + t - p - 1;
    S_2 = beta2 + sum(e2.^2)';      
    Sidraw = gamrnd(S_1./2,2./S_2);
    Sig2 = 1./Sidraw;
     
    
    % Create the VAR covariance matrix H(t). It holds that:
    %           A(t) x H(t) x A(t)' = SIGMA(t) x SIGMA(t) '
    Ht = zeros(M,M,t);
    Htsd = zeros(M,M,t);
    for i = 1:t
        inva = inv(capAt((i-1)*M+1:i*M,:));
        stem = diag(sigt(i,:));
        Hsd = inva*stem;
        Hdraw = Hsd*Hsd';
        Ht(:,:,i) = Hdraw;  % H(t)
        Htsd(:,:,i) = Hsd;  % Cholesky of H(t)
    end