diff --git a/matlab/PosteriorIRF.m b/matlab/PosteriorIRF.m
index 6b5726642..37d196d27 100644
--- a/matlab/PosteriorIRF.m
+++ b/matlab/PosteriorIRF.m
@@ -2,7 +2,6 @@ function PosteriorIRF()
 % stephane.adjemian@ens.fr [09-25-2005]
 global options_ estim_params_ oo_ M_
 
-
 nvx  = estim_params_.nvx;
 nvn  = estim_params_.nvn;
 ncx  = estim_params_.ncx;
@@ -10,7 +9,12 @@ ncn  = estim_params_.ncn;
 np   = estim_params_.np ;
 npar = nvx+nvn+ncx+ncn+np;
 offset = npar-np;
-
+%%
+MaxNumberOfPlotPerFigure = 4;% The square root must be an integer!
+nn = sqrt(MaxNumberOfPlotPerFigure);
+%%
+CheckPath('Plots\IRFs');
+CheckPath('metropolis\IRFs');
 DirectoryName = CheckPath('metropolis');
 load([ DirectoryName '\'  M_.fname '_mh_history'])
 FirstMhFile = record.KeepedDraws.FirstMhFile;
@@ -20,26 +24,28 @@ TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
 NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
 clear record;
 MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
-MAX_nirfs = ceil(options_.MaxNumberOfBytes/(options_.irf*length(oo_.steady_state)*M_.exo_nbr)/8);
-
+MAX_nirfs = ceil(options_.MaxNumberOfBytes/(options_.irf*length(oo_.steady_state)*M_.exo_nbr)/8)+50;
+%%
 B = round(0.25*NumberOfDraws);
-
 if B <= MAX_nirfs
   stock_irf = zeros(options_.irf,size(M_.endo_names,1),M_.exo_nbr,B);
 elseif nvn & B > MAX_nirfs
   stock_irf = zeros(options_.irf,size(M_.endo_names,1),M_.exo_nbr,MAX_nirfs);
 end
-
+%%
 irun = 0;
 ifil = 1;
-
 h = waitbar(0,'Bayesian IRFs...');
-
+if B >= MAX_nirfs 
+  stock_irf = zeros(options_.irf,M_.endo_nbr,M_.exo_nbr,MAX_nirfs);
+else
+  stock_irf = zeros(options_.irf,M_.endo_nbr,M_.exo_nbr,B);
+end
 for b=1:B
   irun = irun+1;
   deep = GetOneDraw(NumberOfDraws,FirstMhFile,LastMhFile,FirstLine,MAX_nruns,DirectoryName);
   M_.params(estim_params_.param_vals(:,1)) = deep(offset+1:end);
-  dr = resol(oo_.steady_state,0)
+  dr = resol(oo_.steady_state,0);
   if nvx
     ip = 1;
     for i=1:nvx
@@ -66,9 +72,9 @@ for b=1:B
       if options_.relative_irf
 	y = 100*y/cs(i,i);
       end
-      for j = 1:size(M_.endo_names,1)
+      for j = 1:M_.endo_nbr%size(M_.endo_names,1)
 	if max(y(j,:)) - min(y(j,:)) > 1e-10 
-	  stock_irf(:,j,i,b) = transpose(y(j,:));
+	  stock_irf(:,j,i,irun) = transpose(y(j,:));
 	end
       end
     end
@@ -85,13 +91,31 @@ for b=1:B
 end
 ifil = ifil-1;
 close(h)
-
+varlist = options_.varlist;
+if isempty(varlist)
+  varlist = M_.endo_names;
+  SelecVariables = transpose(1:M_.endo_nbr);
+  nvar = M_.endo_nbr;
+else
+  nvar = size(varlist,1);
+  SelecVariables = [];
+  for i=1:nvar
+    if ~isempty(strmatch(varlist(i,:),M_.endo_names,'exact'))
+      SelecVariables = [SelecVariables;strmatch(varlist(i,:),M_.endo_names,'exact')];
+    end
+  end
+end
 MeanIRF = zeros(options_.irf,nvar,M_.exo_nbr);
 MedianIRF = zeros(options_.irf,nvar,M_.exo_nbr);
 StdIRF = zeros(options_.irf,nvar,M_.exo_nbr);
 DistribIRF = zeros(options_.irf,9,nvar,M_.exo_nbr);
 HPDIRF = zeros(options_.irf,2,nvar,M_.exo_nbr);
-
+if options_.TeX
+  varlist_TeX = [];
+  for i=1:nvar
+    varlist_TeX = strvcat(varlist_TeX,M_.endo_names_tex(SelecVariables(i),:));
+  end
+end
 fprintf('MH: Posterior IRFs...\n');
 tit(M_.exo_names_orig_ord,:) = M_.exo_names;
 for i = 1:M_.exo_nbr
@@ -106,8 +130,8 @@ for i = 1:M_.exo_nbr
 	StartLine = StartLine+DeProfundis;
       end
       [MeanIRF(k,j,i),MedianIRF(k,j,i),VarIRF(k,j,i),HPDIRF(k,:,j,i),DistribIRF(k,:,j,i)] = posterior_moments(tmp,0);
-      disp(['    Variable: ' deblank(M_.endo_names(SelecVariables(j),:)) ', orthogonalized shock to ' deblank(tit(i,:))])	
     end
+    disp(['    Variable: ' deblank(M_.endo_names(SelecVariables(j),:)) ', orthogonalized shock to ' deblank(tit(i,:))])
   end
 end  
 clear stock_irf;
@@ -122,179 +146,79 @@ for i = 1:M_.exo_nbr
     eval(['oo_.PosteriorIRF.HPDsup.' name ' = HPDIRF(:,2,j,i);']);
   end
 end
-
-
-return
-
 %%
 %% 	Finally i build the plots.
 %%
-
-if TeX
-  fidTeX = fopen([M_.fname '_BayesianIRF.TeX'],'w');
+if options_.TeX
+  fidTeX = fopen([M_.dname '\Plots\IRFs\' M_.fname '_BayesianIRF.TeX'],'w');
   fprintf(fidTeX,'%% TeX eps-loader file generated by PosteriorIRF.m (Dynare).\n');
   fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
   fprintf(fidTeX,' \n');
+  titTeX(M_.exo_names_orig_ord,:) = M_.exo_names_tex;
 end
-
-if TeX; titTeX(M_.exo_name_orig_ord,:) = M_.exo_names_tex; end;
-
-% $$$ for i=1:M_.exo_nbr
-% $$$   number_of_plots_to_draw = 0;
-% $$$   index = [];
-% $$$   for j=1:nvar
-% $$$     if MeanIRF(1,j,i)
-% $$$       number_of_plots_to_draw = number_of_plots_to_draw + 1;
-% $$$       index = cat(1,index,j);
-% $$$     end
-% $$$   end
-% $$$   [nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_plots_to_draw);	
-% $$$   if nbplt == 1
-% $$$     if options_.relative_irf
-% $$$       hh = figure('Name',['Relative response to orthogonalized' ...
-% $$$ 			  ' shock to ' tit(i,:)]);
-% $$$     else
-% $$$       hh = figure('Name',['Orthogonalized shock to ' tit(i, ...
-% $$$ 						  :)]);
-% $$$     end
-% $$$     NAMES = [];
-% $$$     if TeX; TEXNAMES = []; end;
-% $$$     for j=1:number_of_plots_to_draw
-% $$$       set(0,'CurrentFigure',hh)
-% $$$       subplot(nr,nc,j);
-% $$$       plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
-% $$$       hold on
-% $$$       for k = 1:9
-% $$$ 	plot(1:options_.irf,DistribIRF(:,index(j),i,k),'-g','linewidth',0.5)
-% $$$       end
-% $$$       plot(1:options_.irf,MeanIRF(:,index(j),i),'-k','linewidth',1)
-% $$$       xlim([1 options_.irf]);
-% $$$       hold off
-% $$$       name    = deblank(M_.endo_names(SelecVariables(index(j)),:));
-% $$$       NAMES   = strvcat(NAMES,name);
-% $$$       if TeX
-% $$$ 	texname = deblank(M_.endo_names_tex(SelecVariables(index(j)),:));
-% $$$ 	TEXNAMES   = strvcat(TEXNAMES,['$' texname '$']);
-% $$$       end
-% $$$       title(name,'Interpreter','none')
-% $$$     end
-% $$$     eval(['print -depsc2 ' M_.fname '_Bayesian_IRF_' deblank(tit(i,:))]);
-% $$$     eval(['print -dpdf ' M_.fname  '_Bayesian_IRF_' deblank(tit(i,:))]);
-% $$$     saveas(hh,[M_.fname  '_Bayesian_IRF_' deblank(tit(i,:)) '.fig']);
-% $$$     if options_.nograph, close(hh), end
-% $$$     if TeX
-% $$$       fprintf(fidTeX,'\\begin{figure}[H]\n');
-% $$$       for jj = 1:number_of_plots_to_draw
-% $$$ 	fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{%s}\n'],deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:)));
-% $$$       end    
-% $$$       fprintf(fidTeX,'\\centering \n');
-% $$$       fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_Bayesian_IRF_%s}\n',M_.fname,deblank(tit(i,:)));
-% $$$       if options_.relative_irf
-% $$$ 	fprintf(fidTeX,['\\caption{Bayesian relative' ...
-% $$$ 			' IRF.}']);
-% $$$       else
-% $$$ 	fprintf(fidTeX,'\\caption{Bayesian IRF.}');
-% $$$       end
-% $$$       fprintf(fidTeX,'\\label{Fig:BayesianIRF:%s}\n',deblank(tit(i,:)));
-% $$$       fprintf(fidTeX,'\\end{figure}\n');
-% $$$       fprintf(fidTeX,' \n');
-% $$$     end    
-% $$$   elseif nbplt > 1
-% $$$     for fig = 1:nbplt-1
-% $$$       if options_.relative_irf
-% $$$ 	hh = figure('Name',['Relative response to orthogonalized' ...
-% $$$ 			    ' shock to ' tit(i,:) ' figure ' int2str(fig) '.']);
-% $$$       else
-% $$$ 	hh = figure('Name',['Orthogonalized shock to ' tit(i,:) ...
-% $$$ 			    ' figure ' int2str(fig) '.']);
-% $$$       end
-% $$$       NAMES = [];
-% $$$       if TeX; TEXNAMES = []; end;
-% $$$       for j=1:nstar
-% $$$ 	jj = (fig-1)*nstar + j;
-% $$$ 	subplot(nr,nc,j);
-% $$$ 	plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
-% $$$ 	hold on
-% $$$ 	for k = 1:9
-% $$$ 	  plot(1:options_.irf,DistribIRF(:,index(jj),i,k),'-g','linewidth',0.5)
-% $$$ 	end
-% $$$ 	plot(1:options_.irf,MeanIRF(:,index(jj),i),'-k','linewidth',1)
-% $$$ 	xlim([1 options_.irf]);
-% $$$ 	hold off
-% $$$ 	name    = deblank(M_.endo_names(SelecVariables(index(jj)),:));
-% $$$ 	NAMES   = strvcat(NAMES,name);
-% $$$ 	if TeX
-% $$$ 	  texname = deblank(M_.endo_names_tex(SelecVariables(index(jj)),:));
-% $$$ 	  TEXNAMES   = strvcat(TEXNAMES,['$' texname '$']);
-% $$$ 	end
-% $$$ 	title(name,'Interpreter','none')
-% $$$       end
-% $$$       eval(['print -depsc2 ' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) int2str(fig)]);
-% $$$       eval(['print -dpdf ' M_.fname  '_Bayesian_IRF_' deblank(tit(i,:)) int2str(fig)]);
-% $$$       saveas(hh,[M_.fname  '_Bayesian_IRF_' deblank(tit(i,:)) int2str(fig) '.fig']);
-% $$$       if options_.nograph, close(hh), end
-% $$$       if TeX
-% $$$ 	fprintf(fidTeX,'\\begin{figure}[H]\n');
-% $$$ 	for jj = 1:nstar
-% $$$ 	  fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{%s}\n'],deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:)));
-% $$$ 	end    
-% $$$ 	fprintf(fidTeX,'\\centering \n');
-% $$$ 	fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_Bayesian_IRF_%s%s}\n',M_.fname,deblank(tit(i,:)),int2str(fig));
-% $$$ 	if options_.relative_irf == 1
-% $$$ 	  fprintf(fidTeX,['\\caption{Bayesian relative' ...
-% $$$ 			  ' IRF.}']);
-% $$$ 	else
-% $$$ 	  fprintf(fidTeX,'\\caption{Bayesian IRF.}');
-% $$$ 	end
-% $$$ 	fprintf(fidTeX,'\\label{Fig:BayesianIRF:%s:%s}\n',deblank(tit(i,:)), int2str(fig));
-% $$$ 	fprintf(fidTeX,'\\end{figure}\n');
-% $$$ 	fprintf(fidTeX,' \n');
-% $$$       end    
-% $$$     end
-% $$$     hh = figure('Name',['Orthogonalized shock to ' tit(i,:) ' figure ' int2str(nbplt) '.']);
-% $$$     NAMES = [];
-% $$$     if TeX; TEXNAMES = []; end;
-% $$$     for j=1:number_of_plots_to_draw -(nbplt-1)*nstar
-% $$$       jj = (nbplt-1)*nstar + j;
-% $$$       subplot(nr,nc,j);
-% $$$       plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
-% $$$       hold on
-% $$$       for k = 1:9
-% $$$ 	plot(1:options_.irf,DistribIRF(:,index(jj),i,k),'-g','linewidth',0.5)
-% $$$       end
-% $$$       plot(1:options_.irf,MeanIRF(:,index(jj),i),'-k','linewidth',1)
-% $$$       xlim([1 options_.irf]);
-% $$$       hold off
-% $$$       name    = deblank(M_.endo_names(SelecVariables(index(jj)),:));
-% $$$       NAMES   = strvcat(NAMES,name);
-% $$$       if TeX
-% $$$ 	texname = deblank(M_.endo_names_tex(SelecVariables(index(jj)),:));
-% $$$ 	TEXNAMES   = strvcat(TEXNAMES,['$' texname '$']);
-% $$$       end
-% $$$       title(name,'Interpreter','none')
-% $$$     end
-% $$$     eval(['print -depsc2 ' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) int2str(nbplt)]);
-% $$$     eval(['print -dpdf ' M_.fname  '_Bayesian_IRF_' deblank(tit(i,:)) int2str(nbplt)]);
-% $$$     saveas(hh,[M_.fname  '_Bayesian_IRF_' deblank(tit(i,:)) int2str(nbplt) '.fig']);
-% $$$     if options_.nograph, close(hh), end
-% $$$     if TeX
-% $$$       fprintf(fidTeX,'\\begin{figure}[H]\n');
-% $$$       for jj = 1:nstar
-% $$$ 	fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{%s}\n'],deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:)));
-% $$$       end    
-% $$$       fprintf(fidTeX,'\\centering \n');
-% $$$       fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_Bayesian_IRF_%s%s}\n',M_.fname,deblank(tit(i,:)),int2str(nbplt));
-% $$$       fprintf(fidTeX,'\\caption{Bayesian IRF.}');
-% $$$       fprintf(fidTeX,'\\label{Fig:BayesianIRF:%s:%s}\n',deblank(tit(i,:)), int2str(nbplt));
-% $$$       fprintf(fidTeX,'\\end{figure}\n');
-% $$$       fprintf(fidTeX,' \n');
-% $$$     end    
-% $$$   else % nbplt = 0
-% $$$     disp('There''s nothing to plot here!')
-% $$$   end
-% $$$ end
-% $$$ if TeX
-% $$$   fprintf(fidTeX,'%% End of TeX file.\n');
-% $$$   fclose(fidTeX);
-% $$$ end
-% $$$ fprintf('MH: Posterior IRFs, done!\n');
+%%
+figunumber = 0;
+subplotnum = 0;
+for i=1:M_.exo_nbr
+  NAMES = [];
+  if options_.TeX 
+    TEXNAMES = []; 
+  end
+  for j=1:nvar
+    if max(abs(MeanIRF(:,j,i))) > 10^(-6)
+      subplotnum = subplotnum+1;
+      if subplotnum == 1 & options_.relative_irf
+	hh = figure('Name',['Relative response to orthogonalized shock to ' tit(i,:)]);
+      elseif subplotnum == 1 & ~options_.relative_irf
+	hh = figure('Name',['Orthogonalized shock to ' tit(i,:)]);
+      end
+      set(0,'CurrentFigure',hh)
+      subplot(nn,nn,subplotnum);
+      plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
+      hold on
+      for k = 1:9
+ 	plot(1:options_.irf,DistribIRF(:,k,j,i),'-g','linewidth',0.5)
+      end
+      plot(1:options_.irf,MeanIRF(:,j,i),'-k','linewidth',1)
+      xlim([1 options_.irf]);
+      hold off
+      name = deblank(varlist(j,:));
+      NAMES = strvcat(NAMES,name);
+      if options_.TeX
+ 	texname = deblank(varlist_TeX(j,:));
+ 	TEXNAMES = strvcat(TEXNAMES,['$' texname '$']);
+      end
+      title(name,'Interpreter','none')
+    end
+    if subplotnum == MaxNumberOfPlotPerFigure | j == nvar  
+      eval(['print -depsc2 ' M_.dname '\Plots\IRFs\'  M_.fname '_Bayesian_IRF_' deblank(tit(i,:))]);
+      eval(['print -dpdf ' M_.dname '\Plots\IRFs\' M_.fname  '_Bayesian_IRF_' deblank(tit(i,:))]);
+      saveas(hh,[M_.dname '\Plots\IRFs\' M_.fname  '_Bayesian_IRF_' deblank(tit(i,:)) '.fig']);
+      if options_.nograph, close(hh), end
+      if options_.TeX
+	fprintf(fidTeX,'\\begin{figure}[H]\n');
+	for jj = 1:size(TEXNAMES,1)
+	  fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{%s}\n'],deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:)));
+	end    
+	fprintf(fidTeX,'\\centering \n');
+	fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_Bayesian_IRF_%s}\n',M_.fname,deblank(tit(i,:)));
+	if options_.relative_irf
+	  fprintf(fidTeX,['\\caption{Bayesian relative IRF.}']);
+	else
+	  fprintf(fidTeX,'\\caption{Bayesian IRF.}');
+	end
+	fprintf(fidTeX,'\\label{Fig:BayesianIRF:%s}\n',deblank(tit(i,:)));
+	fprintf(fidTeX,'\\end{figure}\n');
+	fprintf(fidTeX,' \n');
+      end
+      subplotnum = 0;
+      figunumber = figunumber+1;
+    end
+  end% loop over selected endo_var
+end% loop over exo_var  
+%%
+if options_.TeX
+  fprintf(fidTeX,'%% End of TeX file.\n');
+  fclose(fidTeX);
+end
+fprintf('MH: Posterior IRFs, done!\n');
\ No newline at end of file
diff --git a/matlab/dynare_estimation.m b/matlab/dynare_estimation.m
index 5f9f6d16a..d62c81c15 100644
--- a/matlab/dynare_estimation.m
+++ b/matlab/dynare_estimation.m
@@ -471,7 +471,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
     ip = nvx+nvn+ncx+ncn+1;
     for i=1:np
       fprintf(fidTeX,'$%s$ & %s & %7.3f & %6.4f & %8.4f & %7.4f \\\\ \n',...
-	      deblank(estim_params_.tex(i,:)),...
+	      M_.param_names_tex(estim_params_.param_vals(i,1),:),...%deblank(estim_params_.tex(i,:)),...
 	      deblank(pnames(bayestopt_.pshape(ip)+1,:)),...
 	      bayestopt_.pmean(ip),...
 	      estim_params_.param_vals(i,6),...