function irf=LP_gibbs_H(dataMatrix,L,...
  REPS,BURN,SKIP,...
    yd,xd,S0b,T0b,shockS,trend,h,A0save)


idx=BURN+1:SKIP:REPS;
fsize=length(idx);
N=cols(dataMatrix);
irf=zeros(fsize,N);
 Y=dataMatrix;
X=[];
for j=1+h-1:L+h-1
    X=[X lagmatrix(Y,j)];
end

tr=(0:rows(Y)-1)';

X=[X ones(rows(Y),1)];
EX=1;
if trend==1
X=[X tr];
EX=EX+1;
elseif trend==2
    X=[X tr tr.^2];
    EX=EX+2;
end

Y=Y(L+h:end,:);
X=X(L+h:end,:);
T=rows(X);
%MCMC
Sigma=eye(N);
gibbs1=1;
gibbs=1;
while gibbs1<=fsize
[gibbs h]
 bdraw=drawvarnewx_2([Y;yd],[X;xd],Sigma,N,L,EX);

 
varcoef=reshape(bdraw,N*L+EX,N);
resid=Y-X*varcoef;
%draw
%scale=resid'*resid+S0b;
scale=getNW(resid,T,h-1,resid'*resid+S0b);
Sigma=iwishrnd(scale,T+T0b);%iwpq(T+T0b,invpd(scale));


if gibbs>BURN
   
if sum(gibbs==idx)>0

%IRF using MR identification
%IRF using MR identification
A0x=squeeze(A0save(gibbs1,:,:));
irf(gibbs1,:)=shockS'*A0x*varcoef(1:N,:);


gibbs1=gibbs1+1;

end
 
end

gibbs=gibbs+1;
end