diff --git a/matlab/stab_map_.m b/matlab/stab_map_.m
index f845e787d..ae4f26ab1 100644
--- a/matlab/stab_map_.m
+++ b/matlab/stab_map_.m
@@ -55,6 +55,10 @@ nshock = estim_params_.nvx;
 nshock = nshock + estim_params_.nvn;
 nshock = nshock + estim_params_.ncx;
 nshock = nshock + estim_params_.ncn;
+number_of_grid_points = 2^9;      % 2^9 = 512 !... Must be a power of two.
+bandwidth = 0;                    % Rule of thumb optimal bandwidth parameter.
+kernel_function = 'gaussian';     % Gaussian kernel for Fast Fourrier Transform approximaton.  
+%kernel_function = 'uniform';     % Gaussian kernel for Fast Fourrier Transform approximaton.  
 
 
 if nargin==0,
@@ -95,11 +99,13 @@ if fload==0 | nargin<2 | isempty(fload),
     options_.irf=0;
     options_.noprint=1;
     for j=1:Nsam,
-        for i=1:estim_params_.np,
-            evalin('base',[bayestopt_.name{i+nshock}, '= ',sprintf('%0.15e',lpmat(j,i)),';'])
-        end
+%         for i=1:estim_params_.np,
+%             evalin('base',[bayestopt_.name{i+nshock}, '= ',sprintf('%0.15e',lpmat(j,i)),';'])
+%         end
+        M_.params(estim_params_.param_vals(:,1)) = lpmat(j,:)';
         %evalin('base','stoch_simul(var_list_);');
         stoch_simul([]);
+        dr_ = oo_.dr;
         %egg(:,j) = sort(eigenvalues_);
         %egg(:,j) = sort(dr_.eigval);
         if isfield(dr_,'eigval'),
@@ -117,9 +123,9 @@ if fload==0 | nargin<2 | isempty(fload),
     % map stable samples
     ix=[1:Nsam];
     for j=1:Nsam,
-        if abs(egg(dr_.npred+1,j))<=options_.qz_criterium; %1+1.e-5;
+        if abs(egg(dr_.npred,j))>=options_.qz_criterium; %(1-(options_.qz_criterium-1)); %1-1.e-5;
             ix(j)=0;
-        elseif abs(egg(dr_.npred,j))>=options_.qz_criterium; %(1-(options_.qz_criterium-1)); %1-1.e-5;
+        elseif (dr_.nboth | dr_.nfwrd) & abs(egg(dr_.npred+1,j))<=options_.qz_criterium; %1+1.e-5;
             ix(j)=0;
         end
     end
@@ -129,8 +135,14 @@ if fload==0 | nargin<2 | isempty(fload),
     ixx=[1:Nsam];
     for j=1:Nsam,
         %if abs(egg(dr_.npred+1,j))>1+1.e-5 & abs(egg(dr_.npred,j))<1-1.e-5;
-        if abs(egg(dr_.npred+1,j))>options_.qz_criterium & abs(egg(dr_.npred,j))<options_.qz_criterium; %(1-(options_.qz_criterium-1));
-            ixx(j)=0;
+        if (dr_.nboth | dr_.nfwrd),
+            if abs(egg(dr_.npred+1,j))>options_.qz_criterium & abs(egg(dr_.npred,j))<options_.qz_criterium; %(1-(options_.qz_criterium-1));
+                ixx(j)=0;
+            end
+        else
+            if abs(egg(dr_.npred,j))<options_.qz_criterium; %(1-(options_.qz_criterium-1));
+                ixx(j)=0;
+            end
         end
     end
     ixx=ixx(find(ixx));   % unstable params
@@ -143,31 +155,42 @@ end
 delete([fname_,'_stab_*.*']);
 delete([fname_,'_stab_SA_*.*']);
 delete([fname_,'_stab_corr_*.*']);
+delete([fname_,'_unstab_corr_*.*']);
 
+if length(ixx)>0 & length(ixx)<Nsam,
 % Blanchard Kahn
 for i=1:ceil(estim_params_.np/12),
     figure,
     for j=1+12*(i-1):min(estim_params_.np,12*i),
         subplot(3,4,j-12*(i-1))
-        hist(lpmat(ix,j),30)
+        optimal_bandwidth = mh_optimal_bandwidth(lpmat(ix,j),length(ix),bandwidth,kernel_function); 
+        [x1,f1] = kernel_density_estimate(lpmat(ix,j),number_of_grid_points,...
+            optimal_bandwidth,kernel_function);
+        plot(x1, f1,':k','linewidth',2)
+        optimal_bandwidth = mh_optimal_bandwidth(lpmat(ixx,j),length(ixx),bandwidth,kernel_function); 
+        [x1,f1] = kernel_density_estimate(lpmat(ixx,j),number_of_grid_points,...
+            optimal_bandwidth,kernel_function);
+        hold on, plot(x1, f1,'k','linewidth',2)
+
+        %hist(lpmat(ix,j),30)
         title(bayestopt_.name{j+nshock})
     end
     saveas(gcf,[fname_,'_stab_',int2str(i)])
 end
 
-% TFP STEP & Blanchard; 
+% Smirnov test for Blanchard; 
 for i=1:ceil(estim_params_.np/12),
     figure,
     for j=1+12*(i-1):min(estim_params_.np,12*i),
         subplot(3,4,j-12*(i-1))
         if ~isempty(ix),
             h=cumplot(lpmat(ix,j));
-            set(h,'color',[1 1 1], 'linestyle',':')
+            set(h,'color',[0 0 0], 'linestyle',':')
         end
         hold on,
         if ~isempty(ixx),
              h=cumplot(lpmat(ixx,j));
-             set(h,'color',[1 1 1])
+             set(h,'color',[0 0 0])
         end
 %         if exist('kstest2')==2 & length(ixx)>0 & length(ixx)<Nsam,
 %             [H,P,KSSTAT] = kstest2(lpmat(ix,j),lpmat(ixx,j));
@@ -181,7 +204,6 @@ for i=1:ceil(estim_params_.np/12),
 end
 
 
-if length(ixx)>0 & length(ixx)<Nsam,
 disp(' ')
 disp(' ')
 disp('Starting bivariate analysis:')
@@ -189,39 +211,16 @@ disp('Starting bivariate analysis:')
 c0=corrcoef(lpmat(ix,:));
 c00=tril(c0,-1);
 
-ifig=0;
-j2=0;
-npar=estim_params_.np;
-for j=1:npar,
-    i2=find(abs(c00(:,j))>alpha2);
-    if length(i2)>0,
-        for jx=1:length(i2),
-            j2=j2+1;
-            if mod(j2,12)==1,
-                ifig=ifig+1;
-                figure('name',['Correlations in the stable sample ', num2str(ifig)]),
-            end
-            subplot(3,4,j2-(ifig-1)*12)
-            %             bar(c0(i2,j)), 
-            %             set(gca,'xticklabel',bayestopt_.name(i2)), 
-            %             set(gca,'xtick',[1:length(i2)])
-            %plot(stock_par(ixx(nfilt+1:end,i),j),stock_par(ixx(nfilt+1:end,i),i2(jx)),'.w')
-            %hold on, 
-            plot(lpmat(ix,j),lpmat(ix,i2(jx)),'.')
-            xlabel(deblank(estim_params_.param_names(j,:)),'interpreter','none'), 
-            ylabel(deblank(estim_params_.param_names(i2(jx),:)),'interpreter','none'), 
-            title(['cc = ',num2str(c0(i2(jx),j))])
-            if (mod(j2,12)==0) & j2>0,
-                saveas(gcf,[fname_,'_stab_corr_',int2str(ifig)])
-            end
-        end
-    end
-    if (j==(npar)) & j2>0,
-        saveas(gcf,[fname_,'_stab_corr_',int2str(ifig)])
-    end
+    stab_map_2(lpmat(ix,:),alpha2, 1);
+    stab_map_2(lpmat(ixx,:),alpha2, 0);
     
-end
-%close all
+else
+    if length(ixx)==0,
+        disp('All parameter values in the prior ranges are stable!')
+    else
+        disp('All parameter values in the prior ranges are unstable!')        
+    end
+
 end
 
 
diff --git a/matlab/stab_map_2.m b/matlab/stab_map_2.m
new file mode 100644
index 000000000..f6c7a5f9d
--- /dev/null
+++ b/matlab/stab_map_2.m
@@ -0,0 +1,61 @@
+function stab_map_2(x,alpha2,istab)
+% function stab_map_2(x,alpha2,istab)
+%
+% Copyright (C) 2005 Marco Ratto
+% THIS PROGRAM WAS WRITTEN FOR MATLAB BY
+% Marco Ratto,
+% Unit of Econometrics and Statistics AF
+% (http://www.jrc.cec.eu.int/uasa/),
+% IPSC, Joint Research Centre
+% The European Commission,
+% TP 361, 21020 ISPRA(VA), ITALY
+% marco.ratto@jrc.it 
+%
+
+global bayestopt_ estim_params_ dr_ options_ ys_ fname_
+
+c0=corrcoef(x);
+c00=tril(c0,-1);
+if istab,
+    fig_nam_=[fname_,'_stab_corr_'];
+else
+    fig_nam_=[fname_,'_unstab_corr_'];
+end
+
+ifig=0;
+j2=0;
+npar=estim_params_.np;
+for j=1:npar,
+    i2=find(abs(c00(:,j))>alpha2);
+    if length(i2)>0,
+        for jx=1:length(i2),
+            j2=j2+1;
+            if mod(j2,12)==1,
+                ifig=ifig+1;
+                if istab
+                    figure('name',['Correlations in the stable sample ', num2str(ifig)]),
+                else
+                    figure('name',['Correlations in the unstable sample ', num2str(ifig)]),
+                end
+            end
+            subplot(3,4,j2-(ifig-1)*12)
+            %             bar(c0(i2,j)), 
+            %             set(gca,'xticklabel',bayestopt_.name(i2)), 
+            %             set(gca,'xtick',[1:length(i2)])
+            %plot(stock_par(ixx(nfilt+1:end,i),j),stock_par(ixx(nfilt+1:end,i),i2(jx)),'.k')
+            %hold on, 
+            plot(x(:,j),x(:,i2(jx)),'.')
+            xlabel(deblank(estim_params_.param_names(j,:)),'interpreter','none'), 
+            ylabel(deblank(estim_params_.param_names(i2(jx),:)),'interpreter','none'), 
+            title(['cc = ',num2str(c0(i2(jx),j))])
+            if (mod(j2,12)==0) & j2>0,
+                saveas(gcf,[fig_nam_,int2str(ifig)])
+            end
+        end
+    end
+    if (j==(npar)) & j2>0,
+        saveas(gcf,[fig_nam_,int2str(ifig)])
+    end
+    
+end
+%close all