function [DECOMP_FEVD,NEWDECOMP_FEVD,impulses_news] = tvpnewsdraws(impulses_unanticipated,VVVmean,nrep,M,t,maxh,nhor,p)

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

for iter = 1:nrep
    for irep = 1:t  
        % 3) Perform variance decomposition based on IRF
        IRF = squeeze(impulses_unanticipated(:,:,irep,iter));
        [DECOMP,TOT] = vardecomp_eric(IRF,nhor,p,M);
        
        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)'*VVVmean(irep,:)';           
        end
        
        if newshock(1,maxh)<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(iter,irep,:,:) = DECOMP(1:M,1:nhor);       
        NEWDECOMP_FEVD(iter,irep,:,:) = NEWDECOMP(:,1:nhor);
        impulses_news(iter,irep,:,:)= newshock; 
    end    
end