function [ irf_mat,irf_mat_se,upper,lower]= ...
    LPIV_SW_CF1(y,x,z,xx,zz,n_lag,nirf,pos,bands,y_dif_in)

ncontrol_lag = n_lag;
N=cols(y);
dnobs=rows(y);
sLevel=abs(norminv((1-bands/100)/2,0,1));

% difference indicators .. 0: y variable in levels; 1: y variable in first differencs
y_dif = y_dif_in;
% Construct matrix of controls
tmp = [y];
xcontrol = ones(dnobs,1);
for i = 1:ncontrol_lag;
    xcontrol = [xcontrol lag(tmp,i)];
end;
xreg = [x xcontrol];
zreg = [z xcontrol];
nr = size(xreg,2);
NR=size(xx,2);
irf_mat = NaN(nirf+1,size(y,2));    % Matrix for storing IRFS
irf_mat_se = NaN(nirf+1,size(y,2)); % IRF Standard errors (from HAC)
irf_mat_smpl = NaN(nirf+1,2);           % First and Last observation in Sample
for iseries = 1:size(y,2);
  for ilead = 0:nirf;
        xxreg = lag(xx,-ilead);
        zzreg = lag(zz,-ilead);
      if y_dif(iseries) == 0;
         yreg = lag(y(:,iseries),-ilead);
      elseif y_dif(iseries) == 1;;
         yreg = y(:,iseries);
         for j = 1:ilead;
              yreg = yreg+lag(y(:,iseries),-j);
         end;
      end;
      tmp = packr([yreg xreg zreg xxreg zzreg ]);
      yp = tmp(:,1);
      xp = tmp(:,2:nr+1);
      zp = tmp(:,nr+2:end-(NR*2));
      xxp=tmp(:,end-(NR*2)+1:end-NR);
      zzp=tmp(:,end-NR+1:end);
      xp=[xp(:,1) xxp xp(:,2:end)];
       zp=[zp(:,1) zzp zp(:,2:end)];
      %cp = tmp(:,end);
      [bhat,vbhat,se_bhat]=hac_tsls(yp,xp,zp,nirf+1,1);
      irf_mat(ilead+1,iseries) = bhat(pos);
      irf_mat_se(ilead+1,iseries) = se_bhat(pos);
%       if iseries == 1;
%        irf_mat_smpl(ilead+1,1) = cp(1);
%        irf_mat_smpl(ilead+1,2) = cp(end);
%       end;
  end;
end;
irf_lp_2_est = irf_mat;
irf_lp_2_se = irf_mat_se;
irf_lp_2_smpl = irf_mat_smpl;
upper=irf_mat+sLevel*irf_mat_se;
lower=irf_mat-sLevel*irf_mat_se;