function [newshocklowalt,newshockhighalt] = tvpnewspercLR(impulses_unanticipated,IRFmean,M,t,maxh,nhor,p)


IRFlow = zeros(M^2,nhor,t);
IRFhigh = zeros(M^2,nhor,t);
newshocklowalt = zeros(M,maxh,t);
newshockhighalt = zeros(M,maxh,t);
% DECOMP_low = zeros(M^2,nhor,t);
% DECOMP_high = zeros(M^2,nhor,t);
% NEWSDECOMP_low = zeros(M,nhor,t);
% NEWSDECOMP_high = zeros(M,nhor,t);
pp=16;
for irep = 1:t
    if mod(irep,round(t./12))==0           
        disp([num2str(100*(irep/t)) '% completed'])
    end    
    IRFFc = squeeze(impulses_unanticipated(:,:,irep,:));
%     IRFmean = squeeze(median(impulses_unanticipated(:,:,irep,:),4));
    [SEDalt, low, high] = impln(IRFFc,M,maxh,pp);
    
    IRFlow(:,:,irep) = IRFmean(:,1:nhor,irep) - SEDalt;
    IRFhigh(:,:,irep) = IRFmean(:,1:nhor,irep) + SEDalt;

    CClow = zeros(M,M,nhor+1);
    CChigh = zeros(M,M,nhor+1);
    for ii = 1:nhor+1; CClow(:,:,ii) = reshape(IRFlow(:,ii),M,M); end;
    for ii = 1:nhor+1; CChigh(:,:,ii) = reshape(IRFhigh(:,ii),M,M); end;
    SurpriseImpLow = CClow(:,1,:);
    options = optimset('MaxFunEval',2000000,'MaxIter',1000000,'TolX',1e-15);
    [thetap, ~, flag] = fminsearch(@LongRunEffect,sign(CClow(1,2:end,end))'.*ones(M-1,1),options,squeeze(CClow(1,:,end)));
    zetap=[0;thetap/norm(thetap)];
    [NewsImpLow] = ComputeIrfOneShock(CClow,zetap); %news shock
    
    SurpriseImpHigh = CChigh(:,1,:);
    options = optimset('MaxFunEval',2000000,'MaxIter',1000000,'TolX',1e-15);
    [thetap, ~, flag] = fminsearch(@LongRunEffect,sign(CChigh(1,2:end,end))'.*ones(M-1,1),options,squeeze(CChigh(1,:,end)));
    zetap=[0;thetap/norm(thetap)];
    [NewsImpHigh] = ComputeIrfOneShock(CChigh,zetap); %news shock
    
    newshocklowalt(:,:,irep) = squeeze(NewsImpLow(:,:,1:nhor)); 
    newshockhighalt(:,:,irep) = squeeze(NewsImpHigh(:,:,1:nhor));
    
    if newshocklowalt(1,maxh,irep)<0
        newshocklowalt(:,:,irep) = -newshocklowalt(:,:,irep);
    end
    
    if newshockhighalt(1,maxh,irep)<0
        newshockhighalt(:,:,irep) = -newshockhighalt(:,:,irep);
    end
    % calculate variance decomposition
    
%     SIRFlow = newshocklowalt.^2;
%     SIRFhigh = newshockhighalt.^2;
%     for i = 1:M           
%         SIRFlow(i,:) = cumsum(SIRFlow(i,:));
%         SIRFhigh(i,:) = cumsum(SIRFhigh(i,:));
%     end
%         
%     NEWDECOMPlow = zeros(M,maxh);
%     NEWDECOMPhigh = zeros(M,maxh);
%     for i = 1:M
%         for j = 1:maxh           
%             NEWDECOMPlow(i,j) = SIRFlow(i,j)/TOTlow(j,i);
%             NEWDECOMPhigh(i,j) = SIRFhigh(i,j)/TOThigh(j,i); 
%         end
%     end
%     DECOMP_low(:,:,irep) = DECOMPlow(:,1:nhor);
%     DECOMP_high(:,:,irep) = DECOMPhigh(:,1:nhor);
%     NEWSDECOMP_low(:,:,irep) = NEWDECOMPlow(:,1:nhor); 
%     NEWSDECOMP_high(:,:,irep) = NEWDECOMPhigh(:,1:nhor); 


end