function [irf,irf_VAR]=...
    LP_gibbs_0_mc(dataMatrix,L,...
  REPS,BURN,SKIP,...
    yd,xd,S0b,T0b,shockS,trend,VAR_horz,maxtrys,A0x)


idx=BURN+1:SKIP:REPS;
fsize=length(idx);
N=cols(dataMatrix);
irf=zeros(fsize,2,N);
A0save=zeros(fsize,N,N);
irf_VAR=zeros(fsize,VAR_horz,N);

Y=dataMatrix;
tr=(0:rows(Y)-1)';
X=[];
for j=1:L
    X=[X lagmatrix(Y,j)];
end
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+1:end,:);
X=X(L+1:end,:);
T=rows(X);

bols=vec(X\Y);
bold0=bols;
%MCMC
Sigma=eye(N);
gibbs1=1;
gibbs=1;
while gibbs1<=fsize

 %bdraw = getcoefxxxx(Y,X,B0,iS0,Sigma,N,L,EX);
 %[bdraw,problem] = getcoef_2(X,B0,iS0,Sigma,N,L,EX,bols,maxtrys);
 [bdraw,problemb]=drawvarnewx_1([Y;yd],[X;xd],Sigma,N,L,bold0,EX,maxtrys);
 bold0=bdraw;
 
varcoef=reshape(bdraw,N*L+EX,N);
resid=Y-X*varcoef;
%draw
scale=resid'*resid+S0b;
Sigma=iwishrnd(scale,T+T0b);%iwpq(T+T0b,invpd(scale));


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




irf(gibbs1,1,:)=shockS'*A0x;
irf(gibbs1,2,:)=shockS'*A0x*varcoef(1:N,:);
A0save(gibbs1,:,:)=A0x;
%VAR irf
irfv=irfsimx(varcoef,N,L,A0x,shockS',VAR_horz,EX);
irf_VAR(gibbs1,:,:)=irfv;

gibbs1=gibbs1+1;

end
 
end

gibbs=gibbs+1;
end