function [DECOMP_FEVD, NEWDECOMP_FEVD, impulses_unanticipated, impulses_news] = imp_dec_BS_tvp(Bt,SIGMAt,t,p,M,K,nhor,maxh)

DECOMP_FEVD =  zeros(t,M,nhor);
NEWDECOMP_FEVD = zeros(t,M,nhor);
impulses_unanticipated = zeros(t,M,nhor);
impulses_news = zeros(t,M,nhor);

% Extract impulse responses       
for irep = 1:t 
    % 1) Obtain estimated VAR coefficients and convert them to matrices
    % of appropriate dimensions for impulse response analysis   
    bbtemp = Bt(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)];
           
    shock = chol(squeeze(SIGMAt(:,:,irep)))';
 
            
    % 2) Estimate IRF based on Cholesky
    IRF = IRFVAR(B,shock,p,nhor);
           
    % 3) Perform variance decomposition based on IRF   
    [DECOMP,TOT] = vardecomp_eric(IRF,nhor,p,M);    
    % variance decomposition: TOT is total forecast error variance of each   
    % variable at any horizon.  DECOMP is the forecast error variance   
    % attributable to each Choleski identified shock
        
    V9 = IRF((M+1:M:(M-1)*M+1),1:maxh)';
      
    WWW = zeros(M-1);
    for j = 1:maxh         
        WWW = WWW + (maxh+1-j)*V9(j,:)'*V9(j,:);
    end
    
    [Q,QQ] = eig(WWW);       % Decompose WWW
    RQ = diag(QQ);           % Main diagonal has eigenvalues
    [maxv,maxi] = max(RQ);   % Maximum is first eigenvalue   
    VVV = [0 Q(:,maxi)'];      % Take the 1st eigenvector (PC estimate)

    newshock = zeros(M,maxh);
    for i = 1:M
        newshock(i,:) = IRF((0:M:(M-1)*M)+i,1:maxh)'*VVV';           
    end
    
    if newshock(1,40)<0
        newshock = -1*newshock;   
        VVV = -VVV;
    end
    
    % calculate variance decomposition
    SIRF = newshock.^2; % This squares every element of the IRF
    for i = 1:M           
        SIRF(i,:) = cumsum(SIRF(i,:));
    end
    
    NEWDECOMP = zeros(M,maxh);
    for i = 1:M
        for j = 1:maxh           
            NEWDECOMP(i,j) = SIRF(i,j)/TOT(j,i);   
        end
    end
    % calculate fracton of TFP due to surprise plus news
    TOT_TFP = DECOMP(1,1:maxh) + NEWDECOMP(1,1:maxh);
           
    DECOMP_FEVD(irep,:,:) = DECOMP(1:M,1:nhor);
    NEWDECOMP_FEVD(irep,:,:) = NEWDECOMP(:,1:nhor);
    impulses_unanticipated(irep,:,:)= IRF(1:M,1:nhor); 
    impulses_news(irep,:,:)= newshock; 
end