function [newshock,IRFmean,VVVmean] = news(BMEAN,SIGMAMEAN,M,K,p,nhor,maxh)

B = [BMEAN' ; eye(M*(p-1)) zeros(M*(p-1),M)];
% st dev matrix for structural VAR         
shock = chol(SIGMAMEAN)';
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 = [0 Q(:,maxi)'];      % Take the 1st eigenvector (PC estimate)

for i = 1:M
    newshock(i,:) = IRF((0:M:(M-1)*M)+i,1:maxh)'*VVVmean';           
end

if newshock(1,maxh)<0       
    newshock = -1*newshock;   
    VVVmean = -VVVmean;
end

IRFmean = IRF(:,1:nhor);



