function irf = get_irflp(shock,b,A0,N,horizon)
K=horizon;
irf=zeros(K+1,N);
shockS1=shock(1,:)';
irf(1,:)=shockS1'*A0';
for h=1:K
    B=b(:,:,h+1);
    shockS1=shock(h+1,:)';
    irf(h+1,:)=shockS1'*A0'*B(1:N,:); 
end
