function irf=LP_gibbs(Y,L,Horizon,T0,m_t,...
    prior_tightness,REPS,BURN,SKIP,...
    pre_sample,shockS,normS,shockno)
idx=BURN+1:SKIP:REPS;
fsize=length(idx);
k=cols(m_t);

ex=1;
if pre_sample
Y0=Y(1:T0,:);
else
Y0=Y;
end
[prior_mean,s_d]=get_prior(Y0,L,Horizon);

%data
if pre_sample
dataMatrix=Y(T0+1:end,:);
m_t_1=m_t(T0+1:end,:);
else
dataMatrix=Y;
m_t_1=m_t;


end
N=cols(dataMatrix);
irf=zeros(fsize,Horizon+1,N);
A0save=zeros(fsize,N,N);



%H=0
Y=dataMatrix;
X=[];
for j=1:L
    X=[X lagmatrix(Y,j)];
end
X=[X ones(rows(Y),1)];
Y=Y(L+1:end,:);
X=X(L+1:end,:);
M=m_t_1(L+1:end,:);
T=rows(X);
B0=vec(prior_mean{1});
S0=diag(s_d{1})*prior_tightness;
S0b=eye(N);
T0b=N+1;
iS0=invpd(S0);
Sigma=eye(N);
gibbs1=1;
gibbs=1;
while gibbs1<=fsize
 
 bdraw = getcoefxxxx(Y,X,B0,iS0,Sigma,N,L,ex);

 
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 using MR identification
if k>1
VAR1.m1=M;
VAR1.res=resid;
VAR1.Sigma=Sigma;
VAR1.ord=[1 2];
VAR1 = doProxySVARxx(VAR1);
A0x=zeros(N,N);
tmp=VAR1.b1(:,shockno,1)';
tmp=tmp./tmp(normS);
A0x(shockno,:)=tmp;
else
VAR1.m1=M;
VAR1.res=resid;
VAR1.Sigma=Sigma;
VAR1 = doProxySVAR_single_xx(VAR1);
A0x=zeros(N,N);
tmp=VAR1.b1';
tmp=tmp./tmp(normS);
A0x(1:k,:)=tmp;
end
irf(gibbs1,1,:)=shockS'*A0x;
A0save(gibbs1,:,:)=A0x;
gibbs1=gibbs1+1;

end
 
end

gibbs=gibbs+1;
end


%%%
for h=1:Horizon
    h
 Y=dataMatrix;
X=[];
for j=1+h-1:L+h-1
    X=[X lagmatrix(Y,j)];
end
X=[X ones(rows(Y),1)];
Y=Y(L+h:end,:);
X=X(L+h:end,:);
T=rows(X);
B0=vec(prior_mean{h+1});
S0=diag(s_d{h+1})*prior_tightness;
iS0=invpd(S0);
Sigma=eye(N);
gibbs1=1;
gibbs=1;
while gibbs1<=fsize
 
 bdraw = getcoefxxxx(Y,X,B0,iS0,Sigma,N,L,ex);

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

scale=getNW(resid,T,h-1,resid'*resid+S0b);
Sigma=iwishrnd(scale,T+T0b);


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

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

end
 
end

gibbs=gibbs+1;
end













end  %%%%%
