diff --git a/matlab/dynare_sensitivity.m b/matlab/dynare_sensitivity.m
index 1d9eeef7f..d85ca1f22 100644
--- a/matlab/dynare_sensitivity.m
+++ b/matlab/dynare_sensitivity.m
@@ -36,6 +36,10 @@ x0=[];
 
 options_gsa = set_default_option(options_gsa,'datafile',[]);
 options_gsa = set_default_option(options_gsa,'rmse',0);
+options_gsa = set_default_option(options_gsa,'useautocorr',0);
+
+options_gsa = set_default_option(options_gsa,'moment_calibration',options_.endogenous_prior_restrictions.moment);
+options_gsa = set_default_option(options_gsa,'irf_calibration',options_.endogenous_prior_restrictions.irf);
 if isfield(options_gsa,'nograph'),
     options_.nograph=options_gsa.nograph;
 end
@@ -100,7 +104,6 @@ if options_gsa.identification,
     options_gsa = set_default_option(options_gsa,'trans_ident',0);
     options_gsa = set_default_option(options_gsa,'load_ident_files',0);
     options_gsa = set_default_option(options_gsa,'ar',1);
-    options_gsa = set_default_option(options_gsa,'useautocorr',0);
     options_.ar = options_gsa.ar;
     if options_gsa.morris==0,
         disp('The option morris = 0 is no longer supported (Type I errors)')
@@ -239,6 +242,10 @@ end
 % redform_map(namendo, namlagendo, namexo, icomp, pprior, ilog, threshold)
 
 options_.opt_gsa = options_gsa;
+if ~isempty(options_gsa.moment_calibration) || ~isempty(options_gsa.irf_calibration),
+    map_calibration(OutputDirectoryName, M_, options_, oo_, estim_params_);
+end
+
 if options_gsa.identification,
     map_ident_(OutputDirectoryName,options_gsa);
 end
diff --git a/matlab/endogenous_prior_restrictions.m b/matlab/endogenous_prior_restrictions.m
index e28a726ae..d5b62104b 100644
--- a/matlab/endogenous_prior_restrictions.m
+++ b/matlab/endogenous_prior_restrictions.m
@@ -1,8 +1,8 @@
-function info = endogenous_prior_restrictions(T,R,Model,DynareOptions,DynareResults);
-% Check for prior (sign) restrictions on irf's 
+function [info, info_irf, info_moment, data_irf, data_moment] = endogenous_prior_restrictions(T,R,Model,DynareOptions,DynareResults);
+% Check for prior (sign) restrictions on irf's and theoretical moments
 %
 % INPUTS
-%    T          [double]     n*n state space matrix 
+%    T          [double]     n*n state space matrix
 %    R          [double]     n*k matrix of shocks
 %    Model      [structure]
 %    DynareOptions [structure]
@@ -11,6 +11,9 @@ function info = endogenous_prior_restrictions(T,R,Model,DynareOptions,DynareResu
 % OUTPUTS
 %    info     [double]  check if prior restrictions are matched by the
 %                       model and related info
+%    info_irf [double] array of test checks for all individual irf restrictions
+%    info_moment [double] array of test checks for all individual moment restrictions
+%    
 
 % Copyright (C) 2013 Dynare Team
 %
@@ -29,45 +32,136 @@ function info = endogenous_prior_restrictions(T,R,Model,DynareOptions,DynareResu
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-endo_prior_restrictions= DynareOptions.endogenous_prior_restrictions;
-info=0;
-if isempty(endo_prior_restrictions.irf),
-    % no restriction to be checked
-    return
-end
-if DynareOptions.order>1,
-    error('The algorithm for prior (sign) restrictions on irf''s is currently restricted to first-order decision rules')
-	return
-end
-infos=[0 0];
-varlist=Model.endo_names(DynareResults.dr.order_var,:);
-varlist=varlist(DynareResults.dr.restrict_var_list,:);
-NT=1;
-for j=1:size(endo_prior_restrictions.irf,1),
-    NT=max(NT,endo_prior_restrictions.irf{j,3});
-end
-for t=1:NT,
-    RR = T^(t-1)*R;
-    for j=1:size(endo_prior_restrictions.irf,1),
-	    if endo_prior_restrictions.irf{j,3}~=t,
-		   continue,
-		end
-        iendo=strmatch(endo_prior_restrictions.irf{j,1},varlist,'exact');
-        iexo=strmatch(endo_prior_restrictions.irf{j,2},Model.exo_names,'exact');
-        if (RR(iendo,iexo)>endo_prior_restrictions.irf{j,4}(1)) && (RR(iendo,iexo)<endo_prior_restrictions.irf{j,4}(2)),
-            infos(j,:)=[0, 0];
-        else
-            if RR(iendo,iexo)<endo_prior_restrictions.irf{j,4}(1),
-                delt = (RR(iendo,iexo)-endo_prior_restrictions.irf{j,4}(1))^2;
-            else
-                delt = (RR(iendo,iexo)-endo_prior_restrictions.irf{j,4}(2))^2;
-            end            
-            infos(j,:)=[49, delt];
+info=[0 0];
+info_irf=[];
+info_moment=[];
+data_irf=[];
+data_moment=[];
+
+endo_prior_restrictions.irf= DynareOptions.endogenous_prior_restrictions.irf;
+endo_prior_restrictions.moment= DynareOptions.endogenous_prior_restrictions.moment;
+
+if ~isempty(endo_prior_restrictions.irf),
+   data_irf=cell(size(endo_prior_restrictions.irf,1),1);
+    
+    if DynareOptions.order>1,
+        error('The algorithm for prior (sign) restrictions on irf''s is currently restricted to first-order decision rules')
+        return
+    end
+    varlist=Model.endo_names(DynareResults.dr.order_var,:);
+    if isempty(T),
+        [T,R,SteadyState,infox,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults);
+    else % check if T and R are given in the restricted form!!!
+        if size(T,1)<size(varlist,1),
+            varlist=varlist(DynareResults.dr.restrict_var_list,:);
         end
     end
+    NT=1;
+    for j=1:size(endo_prior_restrictions.irf,1),
+        NT=max(NT,max(endo_prior_restrictions.irf{j,3}));
+    end
+    info_irf=ones(size(endo_prior_restrictions.irf,1),2);
+    for t=1:NT,
+        if ~DynareOptions.relative_irf,
+            RR = T^(t-1)*R*diag(sqrt(diag(Model.Sigma_e)));
+        else
+            RR = T^(t-1)*R*100;
+        end
+        for j=1:size(endo_prior_restrictions.irf,1),
+            if endo_prior_restrictions.irf{j,3}~=t,
+                continue,
+            end
+            iendo=strmatch(endo_prior_restrictions.irf{j,1},varlist,'exact');
+            iexo=strmatch(endo_prior_restrictions.irf{j,2},Model.exo_names,'exact');
+            data_irf{j}=[data_irf{j}; [t RR(iendo,iexo)]];
+            if (RR(iendo,iexo)>endo_prior_restrictions.irf{j,4}(1)) && (RR(iendo,iexo)<endo_prior_restrictions.irf{j,4}(2)),
+                info_irf(j,:)=info_irf(j,:).*[0, 0];
+            else
+                if RR(iendo,iexo)<endo_prior_restrictions.irf{j,4}(1),
+                    delt = (RR(iendo,iexo)-endo_prior_restrictions.irf{j,4}(1))^2;
+                else
+                    delt = (RR(iendo,iexo)-endo_prior_restrictions.irf{j,4}(2))^2;
+                end
+                info_irf(j,:)=info_irf(j,:).*[49, delt];
+            end
+        end
+    end
+    if any(info_irf),
+        info=[49,sum(info_irf(:,2))];
+    end
+    
 end
-if any(infos),
-    info=[49,sum(infos(:,2))];
+
+if ~isempty(endo_prior_restrictions.moment),
+    
+    if DynareOptions.order>1,
+        error('The algorithm for prior (sign) restrictions on moments is currently restricted to first-order decision rules')
+        return
+    end
+    
+    data_moment=cell(size(endo_prior_restrictions.moment,1),1);
+    var_list_=endo_prior_restrictions.moment{1,1};
+    for  j=1:size(endo_prior_restrictions.moment,1),
+        tmp=endo_prior_restrictions.moment{j,1};
+        if ~ismember(tmp,cellstr(var_list_)),
+            var_list_ = char(var_list_, tmp);
+        end
+        tmp=endo_prior_restrictions.moment{j,2};
+        if ~ismember(tmp,cellstr(var_list_)),
+            var_list_ = char(var_list_, tmp);
+        end
+    end
+    NTmax=0;
+    NTmin=0;
+    for j=1:size(endo_prior_restrictions.moment,1),
+        NTmax=max(NTmax,max(endo_prior_restrictions.moment{j,3}));
+        NTmin=min(NTmin,min(endo_prior_restrictions.moment{j,3}));
+    end
+    info_moment=ones(size(endo_prior_restrictions.moment,1),2);
+    nvar = size(var_list_,1);
+    ivar=zeros(nvar,1);
+    for i=1:nvar
+        i_tmp = strmatch(var_list_(i,:),Model.endo_names,'exact');
+        if isempty(i_tmp)
+            error (['One of the variable specified does not exist']) ;
+        else
+            ivar(i) = i_tmp;
+        end
+    end
+    DynareOptions.ar = max(abs(NTmin),NTmax);
+    [gamma_y,stationary_vars] = th_autocovariances(DynareResults.dr, ivar, Model, DynareOptions,1);
+    for t=NTmin:NTmax,
+        RR = gamma_y{abs(t)+1};     
+        if t==0,
+            RR = RR./(sqrt(diag(RR))*sqrt(diag(RR))')-eye(nvar)+diag(diag(gamma_y{t+1})); % becomes correlation            
+        end
+        for j=1:size(endo_prior_restrictions.moment,1),
+            if endo_prior_restrictions.moment{j,3}~=t,
+                continue,
+            end
+            iendo1 = strmatch(endo_prior_restrictions.moment{j,1},var_list_,'exact');
+            iendo2 = strmatch(endo_prior_restrictions.moment{j,2},var_list_,'exact');
+            if t>0
+                tmp0 = iendo1;
+                iendo1=iendo2;
+                iendo2=tmp0;
+            end
+            data_moment{j}=[data_moment{j}; [t RR(iendo1,iendo2)]];
+            if (RR(iendo1,iendo2)>endo_prior_restrictions.moment{j,4}(1)) && (RR(iendo1,iendo2)<endo_prior_restrictions.moment{j,4}(2)),
+                info_moment(j,:)=info_moment(j,:).*[0, 0];
+            else
+                if RR(iendo1,iendo2)<endo_prior_restrictions.moment{j,4}(1),
+                    delt = (RR(iendo1,iendo2)-endo_prior_restrictions.moment{j,4}(1))^2;
+                else
+                    delt = (RR(iendo1,iendo2)-endo_prior_restrictions.moment{j,4}(2))^2;
+                end
+                info_moment(j,:)=info_moment(j,:).*[49, delt];
+            end
+        end
+    end
+    if any(info_moment),
+        info=[49, info(2) + sum(info_moment(:,2))];
+    end
 end
 return
 
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 39f9659d1..1bb391056 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -609,6 +609,7 @@ options_.risky_steadystate = 0;
 % endogenous prior
 options_.endogenous_prior = 0;
 options_.endogenous_prior_restrictions.irf={};
+options_.endogenous_prior_restrictions.moment={};
 
 % OSR Optimal Simple Rules
 options_.osr.tolf=1e-7;
diff --git a/matlab/gsa/map_calibration.m b/matlab/gsa/map_calibration.m
new file mode 100644
index 000000000..9d8a56e89
--- /dev/null
+++ b/matlab/gsa/map_calibration.m
@@ -0,0 +1,168 @@
+function map_calibration(OutputDirectoryName, Model, DynareOptions, DynareResults, EstimatedParameters)
+
+% Copyright (C) 2014 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+
+fname_ = Model.fname;
+pnames = Model.param_names(EstimatedParameters.param_vals(:,1),:);
+pvalue_ks = DynareOptions.opt_gsa.pvalue_ks;
+indx_irf = [];
+indx_moment = [];
+
+skipline()
+disp('Sensitivity analysis for calibration criteria')
+
+filetoload=[OutputDirectoryName '/' fname_ '_prior'];
+load(filetoload,'lpmat','lpmat0','istable','iunstable','iindeterm','iwrong' ,'infox')
+if ~isempty(lpmat0),
+    lpmatx=lpmat0(istable,:);
+else
+    lpmatx=[];
+end
+[Nsam, np] = size(lpmat);
+nshock = size(lpmat0,2);
+npT = np+nshock;
+
+nbr_irf_restrictions = size(DynareOptions.endogenous_prior_restrictions.irf,1);
+mat_irf=cell(nbr_irf_restrictions,1);
+for ij=1:nbr_irf_restrictions,
+    mat_irf{ij}=NaN(Nsam,length(DynareOptions.endogenous_prior_restrictions.irf{ij,3}));
+end
+
+nbr_moment_restrictions = size(DynareOptions.endogenous_prior_restrictions.moment,1);
+mat_moment=cell(nbr_moment_restrictions,1);
+for ij=1:nbr_moment_restrictions,
+    mat_moment{ij}=NaN(Nsam,length(DynareOptions.endogenous_prior_restrictions.moment{ij,3}));
+end
+
+irestrictions = [1:Nsam];
+for j=1:Nsam,
+    Model = set_all_parameters([lpmat0(j,:) lpmat(j,:)]',EstimatedParameters,Model);
+    [Tt,Rr,SteadyState,info] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict');
+    if info(1)==0,
+        [info, info_irf, info_moment, data_irf, data_moment]=endogenous_prior_restrictions(Tt,Rr,Model,DynareOptions,DynareResults);
+        if ~isempty(info_irf)
+            for ij=1:nbr_irf_restrictions,
+                mat_irf{ij}(j,:)=data_irf{ij}(:,2)';
+            end
+            indx_irf(j,:)=info_irf(:,1);
+        end
+        if ~isempty(info_moment)
+            for ij=1:nbr_moment_restrictions,
+                mat_moment{ij}(j,:)=data_moment{ij}(:,2)';
+            end
+            indx_moment(j,:)=info_moment(:,1);
+        end
+    else
+        irestrictions(j)=0;
+    end
+end
+irestrictions=irestrictions(find(irestrictions));
+xmat=[lpmat0(irestrictions,:) lpmat(irestrictions,:)];
+skipline()
+
+if ~isempty(indx_irf),
+    indx_irf = indx_irf(irestrictions,:);
+    h1=dyn_figure(DynareOptions,'name','Evaluation of irf restrictions');
+    nrow=ceil(sqrt(nbr_irf_restrictions));
+    ncol=nrow;
+    if nrow*(nrow-1)>nbr_irf_restrictions,
+        ncol=nrow-1;
+    end
+    for ij=1:nbr_irf_restrictions,
+        figure(h1),
+        mat_irf{ij}=mat_irf{ij}(irestrictions,:);
+        subplot(nrow,ncol,ij),
+        for ik=1:size(mat_irf{ij},2),
+            cumplot(mat_irf{ij}(:,ik)),
+            hold all,
+        end
+        %     hist(mat_irf{ij}),
+        a=axis;
+        x1val=max(DynareOptions.endogenous_prior_restrictions.irf{ij,4}(1),a(1));
+        x2val=min(DynareOptions.endogenous_prior_restrictions.irf{ij,4}(2),a(2));
+        hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'k');
+        set(hp,'FaceAlpha',[0.5])
+        hold off,
+        leg = num2str(DynareOptions.endogenous_prior_restrictions.irf{ij,3}(1));
+        if size(mat_irf{ij},2)>1,
+            leg = [leg,':' ,num2str(DynareOptions.endogenous_prior_restrictions.irf{ij,3}(end))];
+        end
+        title([DynareOptions.endogenous_prior_restrictions.irf{ij,1},' vs ',DynareOptions.endogenous_prior_restrictions.irf{ij,2}, '(', leg,')'],'interpreter','none'),
+        
+        indx1 = find(indx_irf(:,ij)==0);
+        indx2 = find(indx_irf(:,ij)~=0);
+        atitle=[DynareOptions.endogenous_prior_restrictions.irf{ij,1},' vs ',DynareOptions.endogenous_prior_restrictions.irf{ij,2}, '(', leg,')'];
+        fprintf(['%4.1f%% of the prior support matches IRF ',atitle,' inside [%4.1f, %4.1f]\n'],length(indx1)/length(irestrictions)*100,DynareOptions.endogenous_prior_restrictions.irf{ij,4})
+        aname=['irf_calib_',int2str(ij)];
+        atitle=['IRF Calib: Parameter(s) driving ',DynareOptions.endogenous_prior_restrictions.irf{ij,1},' vs ',DynareOptions.endogenous_prior_restrictions.irf{ij,2}, '(', leg,')'];
+        [proba, dproba] = stab_map_1(xmat, indx1, indx2, aname, 0);
+        indplot=find(proba<pvalue_ks);
+        if ~isempty(indplot)
+            stab_map_1(xmat, indx1, indx2, aname, 1, indplot, OutputDirectoryName,[],atitle);
+        end
+    end
+    dyn_saveas(h1,[OutputDirectoryName,filesep,fname_,'_irf_restrictions'],DynareOptions);
+    skipline()
+end
+
+if ~isempty(indx_moment)
+    indx_moment = indx_moment(irestrictions,:);
+    h2=dyn_figure(DynareOptions,'name','Evaluation of moment restrictions');
+    nrow=ceil(sqrt(nbr_moment_restrictions));
+    ncol=nrow;
+    if nrow*(nrow-1)>nbr_moment_restrictions,
+        ncol=nrow-1;
+    end
+    for ij=1:nbr_moment_restrictions,
+        figure(h2),
+        mat_moment{ij}=mat_moment{ij}(irestrictions,:);
+        subplot(nrow,ncol,ij),
+        for ik=1:size(mat_moment{ij},2),
+            cumplot(mat_moment{ij}(:,ik)),
+            hold all,
+        end
+        %     hist(mat_moment{ij}),
+        a=axis;
+        x1val=max(DynareOptions.endogenous_prior_restrictions.moment{ij,4}(1),a(1));
+        x2val=min(DynareOptions.endogenous_prior_restrictions.moment{ij,4}(2),a(2));
+        hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'k');
+        set(hp,'FaceAlpha',[0.5])
+        hold off,
+        leg = num2str(DynareOptions.endogenous_prior_restrictions.moment{ij,3}(1));
+        if size(mat_moment{ij},2)>1,
+            leg = [leg,':' ,num2str(DynareOptions.endogenous_prior_restrictions.moment{ij,3}(end))];
+        end
+        title([DynareOptions.endogenous_prior_restrictions.moment{ij,1},' vs ',DynareOptions.endogenous_prior_restrictions.moment{ij,2},'(',leg,')'],'interpreter','none'),
+        
+        indx1 = find(indx_moment(:,ij)==0);
+        indx2 = find(indx_moment(:,ij)~=0);
+        atitle=[DynareOptions.endogenous_prior_restrictions.moment{ij,1},' vs ',DynareOptions.endogenous_prior_restrictions.moment{ij,2}, '(', leg,')'];
+        fprintf(['%4.1f%% of the prior support matches MOMENT ',atitle,' inside [%4.1f, %4.1f]\n'],length(indx1)/length(irestrictions)*100,DynareOptions.endogenous_prior_restrictions.moment{ij,4})
+        aname=['moment_calib_',int2str(ij)];
+        atitle=['MOMENT Calib: Parameter(s) driving ',DynareOptions.endogenous_prior_restrictions.moment{ij,1},' vs ',DynareOptions.endogenous_prior_restrictions.moment{ij,2}, '(', leg,')'];
+        [proba, dproba] = stab_map_1(xmat, indx1, indx2, aname, 0);
+        indplot=find(proba<pvalue_ks);
+        if ~isempty(indplot)
+            stab_map_1(xmat, indx1, indx2, aname, 1, indplot, OutputDirectoryName,[],atitle);
+        end
+    end
+    dyn_saveas(h2,[OutputDirectoryName,filesep,fname_,'_moment_restrictions'],DynareOptions);
+    skipline()
+end
+return
+
diff --git a/matlab/gsa/stab_map_.m b/matlab/gsa/stab_map_.m
index c184c20aa..83ae6e10b 100644
--- a/matlab/gsa/stab_map_.m
+++ b/matlab/gsa/stab_map_.m
@@ -253,7 +253,7 @@ if fload==0,
     iindeterm=zeros(1,Nsam);
     iwrong=zeros(1,Nsam);
     for j=1:Nsam,
-        M_.params(estim_params_.param_vals(:,1)) = lpmat(j,:)';
+        M_ = set_all_parameters([lpmat0(j,:) lpmat(j,:)]',estim_params_,M_);
         %try stoch_simul([]);
         try
             [Tt,Rr,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
@@ -496,7 +496,7 @@ delete([OutputDirectoryName,filesep,fname_,'_',auname,'_corr_*.*']);
 delete([OutputDirectoryName,filesep,fname_,'_',aunstname,'_corr_*.*']);
 delete([OutputDirectoryName,filesep,fname_,'_',aindname,'_corr_*.*']);
 
-if length(iunstable)>0 && length(iunstable)<Nsam,
+if length(iunstable)>0 ,
     fprintf(['%4.1f%% of the prior support gives unique saddle-path solution.\n'],length(istable)/Nsam*100)
     fprintf(['%4.1f%% of the prior support gives explosive dynamics.\n'],(length(iunstable)-length(iwrong)-length(iindeterm) )/Nsam*100)
     if ~isempty(iindeterm),
@@ -504,137 +504,137 @@ if length(iunstable)>0 && length(iunstable)<Nsam,
     end
     if ~isempty(iwrong),
         skipline()
-        disp(['For ',num2str(length(iwrong)/Nsam*100,'%1.3f'),'\% of the prior support dynare could not find a solution.'])
+        disp(['For ',num2str(length(iwrong)/Nsam*100,'%4.1f'),'% of the prior support dynare could not find a solution.'])
         skipline()
         if any(infox==1),
-            disp(['    For ',num2str(length(find(infox==1))/Nsam*100,'%1.3f'),'\% The model doesn''t determine the current variables uniquely.'])
+            disp(['    For ',num2str(length(find(infox==1))/Nsam*100,'%4.1f'),'% The model doesn''t determine the current variables uniquely.'])
         end
         if any(infox==2),
-            disp(['    For ',num2str(length(find(infox==2))/Nsam*100,'%1.3f'),'\% MJDGGES returned an error code.'])
+            disp(['    For ',num2str(length(find(infox==2))/Nsam*100,'%4.1f'),'% MJDGGES returned an error code.'])
         end
         if any(infox==6),
-            disp(['    For ',num2str(length(find(infox==6))/Nsam*100,'%1.3f'),'\% The jacobian evaluated at the deterministic steady state is complex.'])
+            disp(['    For ',num2str(length(find(infox==6))/Nsam*100,'%4.1f'),'% The jacobian evaluated at the deterministic steady state is complex.'])
         end
         if any(infox==19),
-            disp(['    For ',num2str(length(find(infox==19))/Nsam*100,'%1.3f'),'\% The steadystate routine thrown an exception (inconsistent deep parameters).'])
+            disp(['    For ',num2str(length(find(infox==19))/Nsam*100,'%4.1f'),'% The steadystate routine thrown an exception (inconsistent deep parameters).'])
         end
         if any(infox==20),
-            disp(['    For ',num2str(length(find(infox==20))/Nsam*100,'%1.3f'),'\% Cannot find the steady state.'])
+            disp(['    For ',num2str(length(find(infox==20))/Nsam*100,'%4.1f'),'% Cannot find the steady state.'])
         end
         if any(infox==21),
-            disp(['    For ',num2str(length(find(infox==21))/Nsam*100,'%1.3f'),'\% The steady state is complex.'])
+            disp(['    For ',num2str(length(find(infox==21))/Nsam*100,'%4.1f'),'% The steady state is complex.'])
         end
         if any(infox==22),
-            disp(['    For ',num2str(length(find(infox==22))/Nsam*100,'%1.3f'),'\% The steady has NaNs.'])
+            disp(['    For ',num2str(length(find(infox==22))/Nsam*100,'%4.1f'),'% The steady has NaNs.'])
         end
         if any(infox==23),
-            disp(['    For ',num2str(length(find(infox==23))/Nsam*100,'%1.3f'),'\% M_.params has been updated in the steadystate routine and has complex valued scalars.'])
+            disp(['    For ',num2str(length(find(infox==23))/Nsam*100,'%4.1f'),'% M_.params has been updated in the steadystate routine and has complex valued scalars.'])
         end
         if any(infox==24),
-            disp(['    For ',num2str(length(find(infox==24))/Nsam*100,'%1.3f'),'\% M_.params has been updated in the steadystate routine and has some NaNs.'])
+            disp(['    For ',num2str(length(find(infox==24))/Nsam*100,'%4.1f'),'% M_.params has been updated in the steadystate routine and has some NaNs.'])
         end
         if any(infox==30),
-            disp(['    For ',num2str(length(find(infox==30))/Nsam*100,'%1.3f'),'\% Ergodic variance can''t be computed.'])
+            disp(['    For ',num2str(length(find(infox==30))/Nsam*100,'%4.1f'),'% 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).'])
+            disp(['    For ',num2str(length(find(infox==49))/Nsam*100,'%4.1f'),'% The model violates one (many) endogenous prior restriction(s).'])
         end
 
     end
     skipline()
-    % Blanchard Kahn
-    [proba, dproba] = stab_map_1(lpmat, istable, iunstable, aname,0);
-%     indstab=find(dproba>ksstat);
-    indstab=find(proba<pvalue_ks);
-    disp('Smirnov statistics in driving acceptable behaviour')
-    for j=1:length(indstab),
-        disp([M_.param_names(estim_params_.param_vals(indstab(j),1),:),'   d-stat = ', num2str(dproba(indstab(j)),'%1.3f'),'   p-value = ', num2str(proba(indstab(j)),'%1.3f')])
-    end
-    skipline()
-    if ~isempty(indstab)
-        stab_map_1(lpmat, istable, iunstable, aname, 1, indstab, OutputDirectoryName,[],atitle);
-    end
-    ixun=iunstable(find(~ismember(iunstable,[iindeterm,iwrong])));
-    if ~isempty(iindeterm),
-        [proba, dproba] = stab_map_1(lpmat, [1:Nsam], iindeterm, aindetname ,0);
-%         indindet=find(dproba>ksstat);
-        indindet=find(proba<pvalue_ks);
-        disp('Smirnov statistics in driving indeterminacy')
-        for j=1:length(indindet),
-            disp([M_.param_names(estim_params_.param_vals(indindet(j),1),:),'   d-stat = ', num2str(dproba(indindet(j)),'%1.3f'),'   p-value = ', num2str(proba(indindet(j)),'%1.3f')])
+    if length(iunstable)<Nsam
+        % Blanchard Kahn
+        [proba, dproba] = stab_map_1(lpmat, istable, iunstable, aname,0);
+        %     indstab=find(dproba>ksstat);
+        indstab=find(proba<pvalue_ks);
+        disp('Smirnov statistics in driving acceptable behaviour')
+        for j=1:length(indstab),
+            disp([M_.param_names(estim_params_.param_vals(indstab(j),1),:),'   d-stat = ', num2str(dproba(indstab(j)),'%1.3f'),'   p-value = ', num2str(proba(indstab(j)),'%1.3f')])
         end
         skipline()
-        if ~isempty(indindet)
-            stab_map_1(lpmat, [1:Nsam], iindeterm, aindetname, 1, indindet, OutputDirectoryName,[],aindettitle);
+        if ~isempty(indstab)
+            stab_map_1(lpmat, istable, iunstable, aname, 1, indstab, OutputDirectoryName,[],atitle);
         end
-    end
-
-    if ~isempty(ixun),
-        [proba, dproba] = stab_map_1(lpmat, [1:Nsam], ixun, aunstablename,0);
-%         indunst=find(dproba>ksstat);
-        indunst=find(proba<pvalue_ks);
-        disp('Smirnov statistics in driving instability')
-        for j=1:length(indunst),
-            disp([M_.param_names(estim_params_.param_vals(indunst(j),1),:),'   d-stat = ', num2str(dproba(indunst(j)),'%1.3f'),'   p-value = ', num2str(proba(indunst(j)),'%1.3f')])
+        ixun=iunstable(find(~ismember(iunstable,[iindeterm,iwrong])));
+        if ~isempty(iindeterm),
+            [proba, dproba] = stab_map_1(lpmat, [1:Nsam], iindeterm, aindetname ,0);
+            %         indindet=find(dproba>ksstat);
+            indindet=find(proba<pvalue_ks);
+            disp('Smirnov statistics in driving indeterminacy')
+            for j=1:length(indindet),
+                disp([M_.param_names(estim_params_.param_vals(indindet(j),1),:),'   d-stat = ', num2str(dproba(indindet(j)),'%1.3f'),'   p-value = ', num2str(proba(indindet(j)),'%1.3f')])
+            end
+            skipline()
+            if ~isempty(indindet)
+                stab_map_1(lpmat, [1:Nsam], iindeterm, aindetname, 1, indindet, OutputDirectoryName,[],aindettitle);
+            end
         end
+        
+        if ~isempty(ixun),
+            [proba, dproba] = stab_map_1(lpmat, [1:Nsam], ixun, aunstablename,0);
+            %         indunst=find(dproba>ksstat);
+            indunst=find(proba<pvalue_ks);
+            disp('Smirnov statistics in driving instability')
+            for j=1:length(indunst),
+                disp([M_.param_names(estim_params_.param_vals(indunst(j),1),:),'   d-stat = ', num2str(dproba(indunst(j)),'%1.3f'),'   p-value = ', num2str(proba(indunst(j)),'%1.3f')])
+            end
+            skipline()
+            if ~isempty(indunst)
+                stab_map_1(lpmat, [1:Nsam], ixun, aunstablename, 1, indunst, OutputDirectoryName,[],aunstabletitle);
+            end
+        end
+        
+        if ~isempty(iwrong),
+            [proba, dproba] = stab_map_1(lpmat, [1:Nsam], iwrong, awronguniname,0);
+            %         indwrong=find(dproba>ksstat);
+            indwrong=find(proba<pvalue_ks);
+            disp('Smirnov statistics in driving no solution')
+            for j=1:length(indwrong),
+                disp([M_.param_names(estim_params_.param_vals(indwrong(j),1),:),'   d-stat = ', num2str(dproba(indwrong(j)),'%1.3f'),'   p-value = ', num2str(proba(indwrong(j)),'%1.3f')])
+            end
+            skipline()
+            if ~isempty(indwrong)
+                stab_map_1(lpmat, [1:Nsam], iwrong, awronguniname, 1, indwrong, OutputDirectoryName,[],awrongunititle);
+            end
+        end
+        
         skipline()
-        if ~isempty(indunst)
-            stab_map_1(lpmat, [1:Nsam], ixun, aunstablename, 1, indunst, OutputDirectoryName,[],aunstabletitle);
+        disp('Starting bivariate analysis:')
+        
+        c0=corrcoef(lpmat(istable,:));
+        c00=tril(c0,-1);
+        
+        if length(istable)>10,
+            stab_map_2(lpmat(istable,:),alpha2, pvalue_corr, asname, OutputDirectoryName,xparam1,astitle);
         end
-    end
-
-    if ~isempty(iwrong),
-        [proba, dproba] = stab_map_1(lpmat, [1:Nsam], iwrong, awronguniname,0);
-%         indwrong=find(dproba>ksstat);
-        indwrong=find(proba<pvalue_ks);
-        disp('Smirnov statistics in driving no solution')
-        for j=1:length(indwrong),
-            disp([M_.param_names(estim_params_.param_vals(indwrong(j),1),:),'   d-stat = ', num2str(dproba(indwrong(j)),'%1.3f'),'   p-value = ', num2str(proba(indwrong(j)),'%1.3f')])
+        if length(iunstable)>10,
+            stab_map_2(lpmat(iunstable,:),alpha2, pvalue_corr, auname, OutputDirectoryName,xparam1,autitle);
         end
-        skipline()
-        if ~isempty(indwrong)
-            stab_map_1(lpmat, [1:Nsam], iwrong, awronguniname, 1, indwrong, OutputDirectoryName,[],awrongunititle);
+        if length(iindeterm)>10,
+            stab_map_2(lpmat(iindeterm,:),alpha2, pvalue_corr, aindname, OutputDirectoryName,xparam1,aindtitle);
         end
-    end
-
-    skipline()
-    disp('Starting bivariate analysis:')
-
-    c0=corrcoef(lpmat(istable,:));
-    c00=tril(c0,-1);
-
-    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
-    if length(iindeterm)>10,
-        stab_map_2(lpmat(iindeterm,:),alpha2, pvalue_corr, aindname, OutputDirectoryName,xparam1,aindtitle);
-    end
-    if length(ixun)>10,
-        stab_map_2(lpmat(ixun,:),alpha2, pvalue_corr, aunstname, OutputDirectoryName,xparam1,aunsttitle);
-    end
-    if length(iwrong)>10,
-        stab_map_2(lpmat(iwrong,:),alpha2, pvalue_corr, awrongname, OutputDirectoryName,xparam1,awrongtitle);
-    end
-
-    x0=0.5.*(bayestopt_.ub(1:nshock)-bayestopt_.lb(1:nshock))+bayestopt_.lb(1:nshock);
-    x0 = [x0; lpmat(istable(1),:)'];
-    if istable(end)~=Nsam
-        M_.params(estim_params_.param_vals(:,1)) = lpmat(istable(1),:)';
-        [oo_.dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
-        %     stoch_simul([]);
-    end
-else
-    if length(iunstable)==0,
-        disp('All parameter values in the specified ranges give unique saddle-path solution!')
+        if length(ixun)>10,
+            stab_map_2(lpmat(ixun,:),alpha2, pvalue_corr, aunstname, OutputDirectoryName,xparam1,aunsttitle);
+        end
+        if length(iwrong)>10,
+            stab_map_2(lpmat(iwrong,:),alpha2, pvalue_corr, awrongname, OutputDirectoryName,xparam1,awrongtitle);
+        end
+        
         x0=0.5.*(bayestopt_.ub(1:nshock)-bayestopt_.lb(1:nshock))+bayestopt_.lb(1:nshock);
         x0 = [x0; lpmat(istable(1),:)'];
+        if istable(end)~=Nsam
+            M_.params(estim_params_.param_vals(:,1)) = lpmat(istable(1),:)';
+            [oo_.dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
+            %     stoch_simul([]);
+        end
     else
         disp('All parameter values in the specified ranges are not acceptable!')
         x0=[];
     end
+else
+    disp('All parameter values in the specified ranges give unique saddle-path solution!')
+        x0=0.5.*(bayestopt_.ub(1:nshock)-bayestopt_.lb(1:nshock))+bayestopt_.lb(1:nshock);
+        x0 = [x0; lpmat(istable(1),:)'];
 
 end
 
diff --git a/tests/gsa/ls2003a.mod b/tests/gsa/ls2003a.mod
new file mode 100644
index 000000000..fcbd9b75c
--- /dev/null
+++ b/tests/gsa/ls2003a.mod
@@ -0,0 +1,97 @@
+var y y_s R pie dq pie_s de A y_obs pie_obs R_obs;
+varexo e_R e_q e_ys e_pies e_A;
+
+parameters psi1 psi2 psi3 rho_R tau alpha rr k rho_q rho_A rho_ys rho_pies;
+
+psi1 = 1.54;
+psi2 = 0.25;
+psi3 = 0.25;
+rho_R = 0.5;
+alpha = 0.3;
+rr = 2.51;
+k = 0.5;
+tau = 0.5;
+rho_q = 0.4;
+rho_A = 0.2;
+rho_ys = 0.9;
+rho_pies = 0.7;
+
+
+model(linear);
+y = y(+1) - (tau +alpha*(2-alpha)*(1-tau))*(R-pie(+1))-alpha*(tau +alpha*(2-alpha)*(1-tau))*dq(+1) + alpha*(2-alpha)*((1-tau)/tau)*(y_s-y_s(+1))-A(+1);
+pie = exp(-rr/400)*pie(+1)+alpha*exp(-rr/400)*dq(+1)-alpha*dq+(k/(tau+alpha*(2-alpha)*(1-tau)))*y+k*alpha*(2-alpha)*(1-tau)/(tau*(tau+alpha*(2-alpha)*(1-tau)))*y_s;
+pie = de+(1-alpha)*dq+pie_s;
+R = rho_R*R(-1)+(1-rho_R)*(psi1*pie+psi2*(y+alpha*(2-alpha)*((1-tau)/tau)*y_s)+psi3*de)+e_R;
+dq = rho_q*dq(-1)+e_q;
+y_s = rho_ys*y_s(-1)+e_ys;
+pie_s = rho_pies*pie_s(-1)+e_pies;
+A = rho_A*A(-1)+e_A;
+y_obs = y-y(-1)+A;
+pie_obs = 4*pie;
+R_obs = 4*R;
+end;
+
+shocks;
+	var e_R = 1.25^2;
+	var e_q = 2.5^2;
+	var e_A = 1.89;
+	var e_ys = 1.89;
+	var e_pies = 1.89;
+end;
+
+varobs y_obs R_obs pie_obs dq de;
+
+estimated_params;
+psi1 , gamma_pdf,1.5,0.5;
+psi2 , gamma_pdf,0.25,0.125;
+psi3 , gamma_pdf,0.25,0.125;
+rho_R ,beta_pdf,0.5,0.2;
+alpha ,beta_pdf,0.3,0.1;
+rr ,gamma_pdf,2.5,1;
+k , gamma_pdf,0.5,0.25;
+tau ,gamma_pdf,0.5,0.2;
+rho_q ,beta_pdf,0.4,0.2;
+rho_A ,beta_pdf,0.5,0.2;
+rho_ys ,beta_pdf,0.8,0.1;
+rho_pies,beta_pdf,0.7,0.15;
+stderr e_R,inv_gamma_pdf,1.2533,0.6551;
+stderr e_q,inv_gamma_pdf,2.5066,1.3103;
+stderr e_A,inv_gamma_pdf,1.2533,0.6551;
+stderr e_ys,inv_gamma_pdf,1.2533,0.6551;
+stderr e_pies,inv_gamma_pdf,1.88,0.9827;
+end;
+
+// endogenous prior restrictions
+options_.relative_irf=1;
+irf_calibration;
+y(1:4), e_ys, [ -50 50]; //[first year response]
+//y(1:4), e_ys, [-inf -50]; //[first year response]
+@#for ilag in 21:40
+R_obs(@{ilag}), e_ys, [0 6]; //[response after 4th year to 10th year]
+@#endfor
+end;
+
+/*
+irf_calibration;
+y(1:4), e_ys, [-inf -0.4]; //[first year response]
+@#for ilag in 21:40
+R_obs(@{ilag}), e_ys, [0 0.25]; //[response after 4th year to 10th year]
+@#endfor
+end;
+*/
+
+moment_calibration;
+//y_obs,y_obs, [0.8 1.1]; //[unconditional variance]
+y_obs,y_obs(-(1:4)), +; //[first year acf]
+@#for ilag in 0:1
+y_obs,R_obs(-@{ilag}), -; //[ccf]
+@#endfor
+end;
+
+dynare_sensitivity(prior_range=0);
+/*
+estimation(datafile='data_ca1.m',first_obs=8,nobs=79,mh_nblocks=2, mode_file = ls2003a_mode,
+  prefilter=1,mh_jscale=0.5,mh_replic=5000, mode_compute=0, mh_drop=0.6, bayesian_irf);
+//stoch_simul(irf=40, order=1, relative_irf) y_obs R_obs pie_obs dq de;
+stoch_simul(irf=40, order=1) y_obs R_obs pie_obs dq de;
+*/