function [newshock,IRFmean,VVVmean] = tvpnews(Bt_mean,SIGMA_mean,M,K,p,t,nhor,maxh)


VVVmean = zeros(t,M);
newshock = zeros(M,maxh,t);
IRFmean = zeros(M^2,nhor+1,t);
for irep = 1:t 
    bbtemp = Bt_mean(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)];
    % st dev matrix for structural VAR       
    shock = chol(squeeze(SIGMA_mean(:,:,irep)))';
    IRF = IRFVAR(B,shock,p,nhor);
    [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   
    VVVmean(irep,:) = [0 Q(:,maxi)'];      % Take the 1st eigenvector (PC estimate)
       
    for i = 1:M
        newshock(i,:,irep) = IRF((0:M:(M-1)*M)+i,1:maxh)'*VVVmean(irep,:)';           
    end
        
    if newshock(1,maxh,irep)<0
        newshock(:,:,irep) = -1*newshock(:,:,irep);   
        VVVmean(irep,:) = -VVVmean(irep,:);
    end
    
    IRFmean(:,:,irep) = IRF;  
end


