diff --git a/matlab/gsa/filt_mc_.m b/matlab/gsa/filt_mc_.m
index b4bc83af5..ed6729362 100644
--- a/matlab/gsa/filt_mc_.m
+++ b/matlab/gsa/filt_mc_.m
@@ -524,17 +524,17 @@ else
     % plot trade-offs
     if ~options_.nograph
     a00=jet(size(vvarvecm,1));
-    for ix=1:ceil(length(nsnam)/5),
-        if options_.opt_gsa.ppost
-            temp_name='RMSE Posterior Tradeoffs: Log Posterior';
+    if options_.opt_gsa.ppost
+        temp_name='RMSE Posterior Tradeoffs:';
+    else
+        if options_.opt_gsa.pprior
+            temp_name='RMSE Prior Tradeoffs:';
         else
-            if options_.opt_gsa.pprior
-                temp_name='RMSE Prior Tradeoffs: Log Posterior';
-            else
-                temp_name='RMSE MC Tradeoffs: Log Posterior';
-            end
-        end        
-        hh = dyn_figure(options_,'name',[temp_name,' ',int2str(ix)]);
+            temp_name='RMSE MC Tradeoffs;';
+        end
+    end
+    for ix=1:ceil(length(nsnam)/5),
+        hh = dyn_figure(options_,'name',[temp_name,' observed variables ',int2str(ix)]);
         for j=1+5*(ix-1):min(size(snam2,1),5*ix),
             subplot(2,3,j-5*(ix-1))
             %h0=cumplot(x(:,nsnam(j)+nshock));
@@ -617,7 +617,7 @@ else
                 fnam = ['rmse_mc_',deblank(vvarvecm(i,:))];
             end
         end
-        stab_map_2(x(ixx(1:nfilt0(i),i),:),alpha2,pvalue,fnam, OutDir);
+        stab_map_2(x(ixx(1:nfilt0(i),i),:),alpha2,pvalue,fnam, OutDir,[],[temp_name ' observed variable ' deblank(vvarvecm(i,:))]);
         
         %     [pc,latent,explained] = pcacov(c0);
         %     %figure, bar([explained cumsum(explained)])
diff --git a/matlab/gsa/map_ident_.m b/matlab/gsa/map_ident_.m
index 591c547c2..867dfcfd6 100644
--- a/matlab/gsa/map_ident_.m
+++ b/matlab/gsa/map_ident_.m
@@ -209,7 +209,7 @@ if opt_gsa.morris==1,
     load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAvdec','vdec','ir_vdec','ic_vdec')
   end
   
-  hh = dyn_figure(options_);
+  hh = dyn_figure(options_,'name','Screening identification: variance decomposition');
 %   boxplot(SAvdec,'whis',10,'symbol','r.')
   myboxplot(SAvdec,[],'.',[],10)
   set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
@@ -221,7 +221,7 @@ if opt_gsa.morris==1,
     text(ip,-2,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
   end
   xlabel(' ')
-  title('All variance decomposition')
+  title('Elementary effects variance decomposition')
   dyn_saveas(hh,[OutputDirectoryName,'/',fname_,'_morris_vdec'],options_);
   else
   save([OutputDirectoryName,'/',fname_,'_morris_IDE'],'vdec')
@@ -314,7 +314,7 @@ if opt_gsa.morris==1,
     load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'ac','ir_ac','ic_ac')
   end
   
-  hh=dyn_figure(options_);
+  hh=dyn_figure(options_,'name','Screening identification: theoretical moments');
 %   boxplot(SAcc,'whis',10,'symbol','r.')
   myboxplot(SAcc,[],'.',[],10)
   set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
@@ -326,7 +326,7 @@ if opt_gsa.morris==1,
     text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
   end
   xlabel(' ')
-  title('EET All moments')
+  title('Elementary effects in the moments')
   dyn_saveas(hh,[OutputDirectoryName,'/',fname_,'_morris_moments'],options_);
 %   close(gcf),
 
@@ -709,7 +709,7 @@ if opt_gsa.morris==1,
   else
     load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAnorm','SAmunorm','SAsignorm')
   end
-  hh=dyn_figure(options_); %bar(SAnorm(:,irel))
+  hh=dyn_figure(options_,'name','Screening identification: model'); %bar(SAnorm(:,irel))
 %   boxplot(SAnorm','whis',10,'symbol','r.')
   myboxplot(SAnorm',[],'.',[],10)
   set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
@@ -725,35 +725,35 @@ if opt_gsa.morris==1,
   title('Elementary effects in the model')
   dyn_saveas(hh,[OutputDirectoryName,'/',fname_,'_morris_par'],options_);
 
-  hh=dyn_figure(options_); %bar(SAmunorm(:,irel))
-%   boxplot(SAmunorm','whis',10,'symbol','r.')
-  myboxplot(SAmunorm',[],'.',[],10)
-  set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
-  set(gca,'xlim',[0.5 npT+0.5])
-  set(gca,'ylim',[-1 1])
-  set(gca,'position',[0.13 0.2 0.775 0.7])
-  xlabel(' ')
-  for ip=1:npT,
-    text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
-  end
-  xlabel(' ')
-  title('\mu in the model')
-  dyn_saveas(hh,[OutputDirectoryName,'/',fname_,'_morrismu_par'],options_);
-
-  hh=dyn_figure(options_); %bar(SAsignorm(:,irel))
-%   boxplot(SAsignorm','whis',10,'symbol','r.')
-  myboxplot(SAsignorm',[],'.',[],10)
-  set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
-  set(gca,'xlim',[0.5 npT+0.5])
-  set(gca,'ylim',[0 1])
-  set(gca,'position',[0.13 0.2 0.775 0.7])
-  xlabel(' ')
-  for ip=1:npT,
-    text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
-  end
-  xlabel(' ')
-  title('\sigma in the model')
-  dyn_saveas(hh,[OutputDirectoryName,'/',fname_,'_morrissig_par'],options_);
+%   hh=dyn_figure(options_); %bar(SAmunorm(:,irel))
+% %   boxplot(SAmunorm','whis',10,'symbol','r.')
+%   myboxplot(SAmunorm',[],'.',[],10)
+%   set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
+%   set(gca,'xlim',[0.5 npT+0.5])
+%   set(gca,'ylim',[-1 1])
+%   set(gca,'position',[0.13 0.2 0.775 0.7])
+%   xlabel(' ')
+%   for ip=1:npT,
+%     text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+%   end
+%   xlabel(' ')
+%   title('\mu in the model')
+%   dyn_saveas(hh,[OutputDirectoryName,'/',fname_,'_morrismu_par'],options_);
+% 
+%   hh=dyn_figure(options_); %bar(SAsignorm(:,irel))
+% %   boxplot(SAsignorm','whis',10,'symbol','r.')
+%   myboxplot(SAsignorm',[],'.',[],10)
+%   set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
+%   set(gca,'xlim',[0.5 npT+0.5])
+%   set(gca,'ylim',[0 1])
+%   set(gca,'position',[0.13 0.2 0.775 0.7])
+%   xlabel(' ')
+%   for ip=1:npT,
+%     text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
+%   end
+%   xlabel(' ')
+%   title('\sigma in the model')
+%   dyn_saveas(hh,[OutputDirectoryName,'/',fname_,'_morrissig_par'],options_);
 
   %     figure, bar(SAnorm(:,irel)')
   %     set(gca,'xtick',[1:j0])
diff --git a/matlab/gsa/mc_moments.m b/matlab/gsa/mc_moments.m
index 31fd3bce6..5c2cd034a 100644
--- a/matlab/gsa/mc_moments.m
+++ b/matlab/gsa/mc_moments.m
@@ -17,19 +17,23 @@ function [vdec, cc, ac] = mc_moments(mm, ss, dr)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-global options_ M_
+global options_ M_ estim_params_ oo_
 
   [nr1, nc1, nsam] = size(mm);
+  nobs=size(options_.varobs,1);
   disp('Computing theoretical moments ...')
   h = dyn_waitbar(0,'Theoretical moments ...');
+  vdec = zeros(nobs,M_.exo_nbr,nsam);
+  cc = zeros(nobs,nobs,nsam);
+  ac = zeros(nobs,nobs*options_.ar,nsam);
   
   for j=1:nsam,
-    dr.ghx = mm(:, [1:(nc1-M_.exo_nbr)],j);
-    dr.ghu = mm(:, [(nc1-M_.exo_nbr+1):end], j);
+    oo_.dr.ghx = mm(:, [1:(nc1-M_.exo_nbr)],j);
+    oo_.dr.ghu = mm(:, [(nc1-M_.exo_nbr+1):end], j);
     if ~isempty(ss),
       set_shocks_param(ss(j,:));
     end
-    [vdec(:,:,j), corr, autocorr, z, zz] = th_moments(dr,options_.varobs);
+    [vdec(:,:,j), corr, autocorr, z, zz] = th_moments(oo_.dr,options_.varobs);
     cc(:,:,j)=triu(corr);
     dum=[];
     for i=1:options_.ar
diff --git a/matlab/gsa/redform_map.m b/matlab/gsa/redform_map.m
index 94c56a41d..2a2e9b027 100644
--- a/matlab/gsa/redform_map.m
+++ b/matlab/gsa/redform_map.m
@@ -132,7 +132,7 @@ for j=1:size(anamendo,1)
                 iplo=iplo+1;
                 js=js+1;
                 xdir0 = [adir,filesep,namendo,'_vs_', namexo];
-                if ilog==0,
+                if ilog==0 || ~isempty(threshold),
                     if isempty(threshold)
                         if isempty(dir(xdir0))
                             mkdir(xdir0)
@@ -233,7 +233,7 @@ for j=1:size(anamendo,1)
                 iplo=iplo+1;
                 js=js+1;
                 xdir0 = [adir,filesep,namendo,'_vs_', namlagendo];
-                if ilog==0,
+                if ilog==0 || ~isempty(threshold),
                     if isempty(threshold)
                         if isempty(dir(xdir0))
                             mkdir(xdir0)
diff --git a/matlab/gsa/stab_map_.m b/matlab/gsa/stab_map_.m
index f2d47c78e..c184c20aa 100644
--- a/matlab/gsa/stab_map_.m
+++ b/matlab/gsa/stab_map_.m
@@ -257,6 +257,9 @@ if fload==0,
         %try stoch_simul([]);
         try
             [Tt,Rr,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
+            if info(1)==0,
+                info=endogenous_prior_restrictions(Tt,Rr,M_,options_,oo_);
+            end
             infox(j,1)=info(1);
             if infox(j,1)==0 && ~exist('T'),
                 dr_=oo_.dr;
@@ -468,11 +471,11 @@ if pprior
     aunstablename=[aname, '_unst'];  aunstabletitle='Prior StabMap: Parameter driving explosiveness of solution';
     awronguniname=[aname, '_wrong']; awrongunititle='Prior StabMap: Parameter driving inability to find solution';
     % bivariate
-    auname='prior_unacceptable'; autitle='Prior non-existence of unique stable solution (Unacceptable)';
-    aunstname='prior_unstable'; aunsttitle='Prior explosiveness of solution';
-    aindname='prior_indeterm'; aindtitle='Prior Indeterminacy';
-    awrongname='prior_wrong'; awrongtitle='Prior inability to find solution';
-    asname='prior_stable'; astitle='Prior unique Stable Saddle-Path';
+    auname='prior_unacceptable'; autitle='Prior StabMap: non-existence of unique stable solution (Unacceptable)';
+    aunstname='prior_unstable'; aunsttitle='Prior StabMap: explosiveness of solution';
+    aindname='prior_indeterm'; aindtitle='Prior StabMap: Indeterminacy';
+    awrongname='prior_wrong'; awrongtitle='Prior StabMap: inability to find solution';
+    asname='prior_stable'; astitle='Prior StabMap: unique Stable Saddle-Path';
 else
     % univariate
     aname='mc_stab'; atitle='MC (around posterior mode) StabMap: Parameter driving non-existence of unique stable solution (Unacceptable)';
@@ -480,11 +483,11 @@ else
     aunstablename=[aname, '_unst'];  aunstabletitle='MC (around posterior mode) StabMap: Parameter driving explosiveness of solution';
     awronguniname=[aname, '_wrong']; awrongunititle='MC (around posterior mode) StabMap: Parameter driving inability to find solution';
     % bivariate
-    auname='mc_unacceptable'; autitle='MC (around posterior mode) non-existence of unique stable solution (Unacceptable)';
-    aunstname='mc_unstable'; aunsttitle='MC (around posterior mode) explosiveness of solution';
-    aindname='mc_indeterm';  aindtitle='MC (around posterior mode) Indeterminacy';
-    awrongname='mc_wrong'; awrongtitle='MC (around posterior mode) inability to find solution';
-    asname='mc_stable'; astitle='MC (around posterior mode) Unique Stable Saddle-Path';
+    auname='mc_unacceptable'; autitle='MC (around posterior mode) StabMap: non-existence of unique stable solution (Unacceptable)';
+    aunstname='mc_unstable'; aunsttitle='MC (around posterior mode) StabMap: explosiveness of solution';
+    aindname='mc_indeterm';  aindtitle='MC (around posterior mode) StabMap: Indeterminacy';
+    awrongname='mc_wrong'; awrongtitle='MC (around posterior mode) StabMap: inability to find solution';
+    asname='mc_stable'; astitle='MC (around posterior mode) StabMap: Unique Stable Saddle-Path';
 end
 delete([OutputDirectoryName,filesep,fname_,'_',aname,'_*.*']);
 %delete([OutputDirectoryName,filesep,fname_,'_',aname,'_SA_*.*']);
@@ -533,6 +536,9 @@ if length(iunstable)>0 && length(iunstable)<Nsam,
         if any(infox==30),
             disp(['    For ',num2str(length(find(infox==30))/Nsam*100,'%1.3f'),'\% Ergodic variance can''t be computed.'])
         end
+        if any(infox==49),
+            disp(['    For ',num2str(length(find(infox==49))/Nsam*100,'%1.3f'),'\% The model violates one (many) endogenous prior restriction(s).'])
+        end
 
     end
     skipline()
@@ -597,7 +603,9 @@ if length(iunstable)>0 && length(iunstable)<Nsam,
     c0=corrcoef(lpmat(istable,:));
     c00=tril(c0,-1);
 
-    stab_map_2(lpmat(istable,:),alpha2, pvalue_corr, asname, OutputDirectoryName,xparam1,astitle);
+    if length(istable)>10,
+        stab_map_2(lpmat(istable,:),alpha2, pvalue_corr, asname, OutputDirectoryName,xparam1,astitle);
+    end
     if length(iunstable)>10,
         stab_map_2(lpmat(iunstable,:),alpha2, pvalue_corr, auname, OutputDirectoryName,xparam1,autitle);
     end
diff --git a/matlab/gsa/stab_map_2.m b/matlab/gsa/stab_map_2.m
index e997af919..d4e94f02e 100644
--- a/matlab/gsa/stab_map_2.m
+++ b/matlab/gsa/stab_map_2.m
@@ -86,14 +86,9 @@ for j=1:npar,
                 fprintf(1,'%20s: corrcoef = %7.3f\n',tmp_name,c0(i2(jx),j));
                     
                 if ~options_.nograph,
-                if strcmp(fnam(1:2),'mc')
-                   type='MC (around posterior mode) StabMap: '; 
-                elseif strcmp(fnam(1:5),'prior')
-                   type='Prior StabMap: ';
-                end
                 if mod(j2,12)==1,
                     ifig=ifig+1;
-                    hh=dyn_figure(options_,'name',[type,'Correlations in the ',figtitle,' sample ', num2str(ifig)]);
+                    hh=dyn_figure(options_,'name',[figtitle,' sample bivariate projection ', num2str(ifig)]);
                 end
                 subplot(3,4,j2-(ifig-1)*12)
                 %             bar(c0(i2,j)),
diff --git a/matlab/gsa/th_moments.m b/matlab/gsa/th_moments.m
index 2837fb696..9dbe7f8f6 100644
--- a/matlab/gsa/th_moments.m
+++ b/matlab/gsa/th_moments.m
@@ -35,11 +35,12 @@ function [vdec, corr, autocorr, z, zz] = th_moments(dr,var_list)
     end
   end
   
-  [gamma_y,ivar] = th_autocovariances(dr,ivar,M_, options_);
-  m = dr.ys(ivar);
+  [gamma_y,stationary_vars] = th_autocovariances(dr,ivar,M_, options_);
+  m = dr.ys(ivar(stationary_vars));
 
   
-  i1 = find(abs(diag(gamma_y{1})) > 1e-12);
+%   i1 = find(abs(diag(gamma_y{1})) > 1e-12);
+  i1 = [1:length(ivar)];
   s2 = diag(gamma_y{1});
   sd = sqrt(s2);