function [newshock,IRFmean] = tvpnewsLR(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);
 
                   
    CC = zeros(M,M,nhor+1);
    for ii = 1:nhor+1; CC(:,:,ii) = reshape(IRF(:,ii),M,M); end;                
    SurpriseImp = CC(:,1,:);
    options = optimset('MaxFunEval',2000000,'MaxIter',1000000,'TolX',1e-15);
    [thetap, ~, flag] = fminsearch(@LongRunEffect,sign(CC(1,2:end,end))'.*ones(M-1,1),options,squeeze(CC(1,:,end)));
    zetap=[0;thetap/norm(thetap)];
    [NewsImp] = ComputeIrfOneShock(CC,zetap); %news shock
    
    newshock(:,:,irep) = squeeze(NewsImp(:,:,1:nhor));
    
    IRFmean(:,:,irep) = IRF;  
end


