function [irf,A0save,irf_VAR,shock_save,shocksm,varcoefm]=LP_gibbs_0(dataMatrix,L,m_t,...
  REPS,BURN,SKIP,...
    yd,xd,S0b,T0b,shockS,normS,shockno,trend,ord,VAR_horz,maxtrys)


idx=BURN+1:SKIP:REPS;
fsize=length(idx);
k=cols(m_t);
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,:);
M=m_t(L+1:end,:);
T=rows(X);
shock_save=zeros(fsize,T);

bols=vec(X\Y);
bold0=bols;
%storage
bsave=zeros(fsize,N*L+EX,N);
ssave=zeros(fsize,N,N);
%MCMC
Sigma=eye(N);
gibbs1=1;
gibbs=1;
while gibbs1<=fsize
gibbs
 %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
bsave(gibbs1,:,:)=varcoef;
ssave(gibbs1,:,:)=Sigma;
%IRF using MR identification
tmp=packr([ M resid]);
proxy=tmp(:,1:k);
innovations=tmp(:,k+1:end);

if k>1
% VAR1.m1=M;
% VAR1.res=resid;
VAR1.resx=resid;
VAR1.m1=proxy;
VAR1.res=innovations;
VAR1.Sigma=Sigma;
VAR1.ord=ord;
if k==2
VAR1 = doProxySVARx(VAR1);
else
  VAR1 = doProxySVARx_3(VAR1);
end
A0x=zeros(N,N);
tmp=VAR1.b1(:,shockno,1)';
tmp=tmp./tmp(normS);
A0x(shockno,:)=tmp;
shocks=VAR1.shocks(:,shockno);
else
VAR1.k=k;
VAR1.resx=resid;
VAR1.m1=proxy;
VAR1.res=innovations;
VAR1.Sigma=Sigma;
VAR1 = doProxySVAR_single_xx(VAR1);
A0x=zeros(N,N);
tmp=VAR1.b1';
tmp=tmp./tmp(normS);
A0x(1:k,:)=tmp;
shocks=VAR1.et;
end


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;
shock_save(gibbs1,:)=shocks';

gibbs1=gibbs1+1;

end
 
end

gibbs=gibbs+1;
end

%posterior mean
varcoefm=squeeze(mean(bsave));
sigmam=squeeze(mean(ssave));
residm=Y-X*varcoefm;
tmp=packr([ M residm]);
proxy=tmp(:,1:k);
innovations=tmp(:,k+1:end);

if k>1
VAR1.rex=residm;
VAR1.m1=proxy;
VAR1.res=innovations;
VAR1.Sigma=sigmam;
VAR1.ord=ord;
if k==2
VAR1 = doProxySVARx(VAR1);
else

  VAR1 = doProxySVARx_3(VAR1);
end
shocksm=VAR1.shocks(:,shockno);
else
VAR1.k=k;
VAR1.resx=residm;
VAR1.m1=proxy;
VAR1.res=innovations;
VAR1.Sigma=sigmam;
VAR1 = doProxySVAR_single_xx(VAR1);

shocksm=VAR1.et;
end
