diff --git a/matlab/GetOneDraw.m b/matlab/GetOneDraw.m
index 2d7564cc4a485a2d8d1983e5beea7278a2b35942..e7ea3e7f7785d53a3f659fb116f8d14f73879766 100644
--- a/matlab/GetOneDraw.m
+++ b/matlab/GetOneDraw.m
@@ -36,4 +36,5 @@ switch type
     [xparams, logpost] = metropolis_draw(0);
   case 'prior'
     xparams = prior_draw(0);
+    logpost = evaluate_posterior_kernel(xparams');
 end
\ No newline at end of file
diff --git a/matlab/PosteriorIRF.m b/matlab/PosteriorIRF.m
index b4abf2d099681129da9ecc3f3c9ce787af531a22..691bac6dff1a5ae69466fb5945e8dbd93f47f310 100644
--- a/matlab/PosteriorIRF.m
+++ b/matlab/PosteriorIRF.m
@@ -80,7 +80,11 @@ DirectoryName = CheckPath('Output');
 if strcmpi(type,'posterior')
     MhDirectoryName = CheckPath('metropolis');
 elseif strcmpi(type,'gsa')
-    MhDirectoryName = CheckPath('GSA');
+    if options_.opt_gsa.pprior
+        MhDirectoryName = CheckPath(['GSA' filesep 'prior']);
+    else
+        MhDirectoryName = CheckPath(['GSA' filesep 'mc']);
+    end
 else
     MhDirectoryName = CheckPath('prior');
 end
@@ -91,7 +95,12 @@ if strcmpi(type,'posterior')
     TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
     NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
 elseif strcmpi(type,'gsa')
-    load([ MhDirectoryName filesep  M_.fname '_prior.mat'],'lpmat0','lpmat','istable')
+    RootDirectoryName = CheckPath('gsa');
+    if options_.opt_gsa.pprior
+        load([ RootDirectoryName filesep  M_.fname '_prior.mat'],'lpmat0','lpmat','istable')
+    else
+        load([ RootDirectoryName filesep  M_.fname '_mc.mat'],'lpmat0','lpmat','istable')
+    end
     x=[lpmat0(istable,:) lpmat(istable,:)];
     clear lpmat istable
     NumberOfDraws=size(x,1);
@@ -193,6 +202,7 @@ localVars.MAX_nruns=MAX_nruns;
 localVars.NumberOfIRFfiles_dsge=NumberOfIRFfiles_dsge;
 localVars.NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar;
 localVars.ifil2=ifil2;
+localVars.MhDirectoryName=MhDirectoryName;
 
 % Like sequential execution!
 if isnumeric(options_.parallel),
@@ -238,7 +248,7 @@ if nosaddle
     disp(['PosteriorIRF :: Percentage of discarded posterior draws = ' num2str(nosaddle/(B+nosaddle))]) 
 end
 
-ReshapeMatFiles('irf_dsge')
+ReshapeMatFiles('irf_dsge',type)
 if MAX_nirfs_dsgevar
     ReshapeMatFiles('irf_bvardsge')
 end
diff --git a/matlab/PosteriorIRF_core1.m b/matlab/PosteriorIRF_core1.m
index a4ed0f633d55b83710afd64a8a2c7190b7998323..c0d7e63f8ab33a4b1a5363946b56026f0115a07e 100644
--- a/matlab/PosteriorIRF_core1.m
+++ b/matlab/PosteriorIRF_core1.m
@@ -82,7 +82,7 @@ if whoiam
 end
 
 
-MhDirectoryName = CheckPath('metropolis');
+MhDirectoryName = myinputs.MhDirectoryName;
 
 RemoteFlag = 0;
 
diff --git a/matlab/ReshapeMatFiles.m b/matlab/ReshapeMatFiles.m
index ad44331c7b7ecf38f949ddfc21145a643da8df0f..933ad8550b97fb802e8ff9af5a7e4d776c3c944e 100644
--- a/matlab/ReshapeMatFiles.m
+++ b/matlab/ReshapeMatFiles.m
@@ -54,8 +54,10 @@ else
             MhDirectoryName = [CheckPath('GSA\SCREEN') filesep ];
         elseif options_.opt_gsa.morris==2,
             MhDirectoryName = [CheckPath('GSA\IDENTIF') filesep ];
+        elseif options_.opt_gsa.pprior
+            MhDirectoryName = [CheckPath(['GSA' filesep 'prior']) filesep ];
         else
-            MhDirectoryName = [CheckPath('GSA') filesep ];
+            MhDirectoryName = [CheckPath(['GSA' filesep 'mc']) filesep ];
         end
     else
         MhDirectoryName = [CheckPath('prior') filesep ];
diff --git a/matlab/dynare_estimation_init.m b/matlab/dynare_estimation_init.m
index 33295707bf3d359561aeb00d4c5ed046692e26aa..4fa2a91b05032deae2ec584847402e9133534b14 100644
--- a/matlab/dynare_estimation_init.m
+++ b/matlab/dynare_estimation_init.m
@@ -1,4 +1,4 @@
-function [data,rawdata,xparam1]=dynare_estimation_init(var_list_, dname, gsa_flag)
+function [data,rawdata,xparam1,data_info]=dynare_estimation_init(var_list_, dname, gsa_flag)
 
 % function dynare_estimation_init(var_list_, gsa_flag)
 % preforms initialization tasks before estimation or
@@ -39,6 +39,14 @@ global M_ options_ oo_ estim_params_ bayestopt_
 
 if nargin < 3 || isempty(gsa_flag)
     gsa_flag = 0;
+else
+    %% Decide if a DSGE or DSGE-VAR has to be estimated.
+    if ~isempty(strmatch('dsge_prior_weight',M_.param_names))
+        options_.dsge_var = 1;
+    end
+    
+    var_list_ = check_list_of_variables(options_, M_, var_list_);
+    options_.varlist = var_list_;
 end
 
 options_.lgyidx2varobs = zeros(size(M_.endo_names,1),1);
@@ -262,6 +270,7 @@ if isempty(options_.datafile)
     if gsa_flag
         data = [];
         rawdata = [];
+        data_info = [];
         return
     else
         error('datafile option is missing')
@@ -308,3 +317,44 @@ end
 % Transpose the dataset array.
 data = transpose(rawdata);
 
+if nargout>3
+    %% Compute the steady state: 
+    if options_.steadystate_flag% if the *_steadystate.m file is provided.
+        [ys,tchek] = feval([M_.fname '_steadystate'],...
+                           [zeros(M_.exo_nbr,1);...
+                            oo_.exo_det_steady_state]);
+        if size(ys,1) < M_.endo_nbr 
+            if length(M_.aux_vars) > 0
+                ys = add_auxiliary_variables_to_steadystate(ys,M_.aux_vars,...
+                                                            M_.fname,...
+                                                            zeros(M_.exo_nbr,1),...
+                                                            oo_.exo_det_steady_state,...
+                                                            M_.params,...
+                                                            options_.bytecode);
+            else
+                error([M_.fname '_steadystate.m doesn''t match the model']);
+            end
+        end
+        oo_.steady_state = ys;
+    else% if the steady state file is not provided.
+        [dd,info] = resol(oo_.steady_state,0);
+        oo_.steady_state = dd.ys; clear('dd');
+    end
+    if all(abs(oo_.steady_state(bayestopt_.mfys))<1e-9)
+        options_.noconstant = 1;
+    else
+        options_.noconstant = 0;
+    end
+    
+    [data_index,number_of_observations,no_more_missing_observations] = describe_missing_data(data,gend,n_varobs);
+    missing_value = ~(number_of_observations == gend*n_varobs);
+    
+    initial_estimation_checks(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations);
+
+    data_info.gend = gend;
+    data_info.data = data;
+    data_info.data_index = data_index;
+    data_info.number_of_observations = number_of_observations;
+    data_info.no_more_missing_observations = no_more_missing_observations;
+    data_info.missing_value = missing_value;
+end
diff --git a/matlab/dynare_sensitivity.m b/matlab/dynare_sensitivity.m
index ea9bd6fc50af91c632fbace9b9ca59a73702b99b..be85312cc406ce0526ebbb157686aa9ba380d72a 100644
--- a/matlab/dynare_sensitivity.m
+++ b/matlab/dynare_sensitivity.m
@@ -24,7 +24,6 @@ function x0=dynare_sensitivity(options_gsa)
 global M_ options_ oo_ bayestopt_ estim_params_
 
 fname_ = M_.fname;
-M_.dname = fname_;
 lgy_ = M_.endo_names;
 x0=[];
 
@@ -57,7 +56,9 @@ if ~isempty(options_gsa.datafile) || isempty(bayestopt_),
         options_.loglinear=options_gsa.loglinear;
     end
     options_.mode_compute = 0;
-    [data,rawdata]=dynare_estimation_init([],1);
+    options_.filtered_vars = 1;
+    options_.plot_priors = 0;
+    [data,rawdata,xparam1,data_info]=dynare_estimation_init([],fname_,1);
     % computes a first linear solution to set up various variables
 end
 
@@ -233,7 +234,11 @@ if options_gsa.rmse,
             a=whos('-file',[OutputDirectoryName,'/',fname_,'_mc'],'logpo2');
         end
         if isempty(a),
-            dynare_MC([],OutputDirectoryName);
+%             dynare_MC([],OutputDirectoryName,data,rawdata,data_info);
+            prior_posterior_statistics('gsa',data,data_info.gend,data_info.data_index,data_info.missing_value);
+            if options_.bayesian_irf
+                PosteriorIRF('gsa');
+            end
             options_gsa.load_rmse=0;
             %   else
             %     if options_gsa.load_rmse==0,
@@ -252,6 +257,7 @@ if options_gsa.rmse,
         end
     end
     clear a;
+%     filt_mc_(OutputDirectoryName,data_info);
     filt_mc_(OutputDirectoryName);
 end
 
diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m
index 64a2ebb2dfb0a3da571bc9616a4857122ab06ea9..3cf6c2cadf9c0b8bbbb2db0aa3cc9b1442d352d3 100644
--- a/matlab/prior_posterior_statistics.m
+++ b/matlab/prior_posterior_statistics.m
@@ -65,23 +65,44 @@ if horizon
 end
 maxlag = M_.maximum_endo_lag;
 %%
-DirectoryName = CheckPath('metropolis');
-load([ DirectoryName '/'  M_.fname '_mh_history'])
-FirstMhFile = record.KeepedDraws.FirstMhFile;
-FirstLine = record.KeepedDraws.FirstLine;
-TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); LastMhFile = TotalNumberOfMhFiles;
-TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
-NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
-clear record;
-if ~isempty(options_.subdraws)
-    B = options_.subdraws;
-    if B > NumberOfDraws
-        B = NumberOfDraws;
+if strcmpi(type,'posterior')
+    DirectoryName = CheckPath('metropolis');
+    load([ DirectoryName '/'  M_.fname '_mh_history'])
+    FirstMhFile = record.KeepedDraws.FirstMhFile;
+    FirstLine = record.KeepedDraws.FirstLine;
+    TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); LastMhFile = TotalNumberOfMhFiles;
+    TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
+    NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
+    clear record;
+    if ~isempty(options_.subdraws)
+        B = options_.subdraws;
+        if B > NumberOfDraws
+            B = NumberOfDraws;
+        end
+    else
+        B = min(1200, round(0.25*NumberOfDraws));
+    end
+elseif strcmpi(type,'gsa')
+    RootDirectoryName = CheckPath('GSA');
+    if options_.opt_gsa.pprior
+        DirectoryName = CheckPath(['GSA',filesep,'prior']);
+        load([ RootDirectoryName filesep  M_.fname '_prior.mat'],'lpmat0','lpmat','istable')
+    else
+        DirectoryName = CheckPath(['GSA',filesep,'mc']);
+        load([ RootDirectoryName filesep  M_.fname '_mc.mat'],'lpmat0','lpmat','istable')
+    end
+    x=[lpmat0(istable,:) lpmat(istable,:)];
+    clear lpmat istable
+    NumberOfDraws=size(x,1);
+    B=NumberOfDraws; 
+elseif strcmpi(type,'prior')
+    DirectoryName = CheckPath('prior');
+    if ~isempty(options_.subdraws)
+        B = options_.subdraws;
+    else
+        B = 1200;
     end
-else
-    B = min(1200, round(0.25*NumberOfDraws));
 end
-
 %%
 MAX_nruns = min(B,ceil(MaxNumberOfBytes/(npar+2)/8));
 MAX_nsmoo = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*gend)/8));
@@ -176,7 +197,7 @@ end
 localVars.MAX_nruns=MAX_nruns;
 localVars.MAX_momentsno = MAX_momentsno;
 localVars.ifil=ifil;
-
+localVars.DirectoryName = DirectoryName;
 
 if strcmpi(type,'posterior'),
     b=0;
@@ -184,11 +205,11 @@ if strcmpi(type,'posterior'),
         b = b + 1;
         [x(b,:), logpost(b,1)] = GetOneDraw(type);
     end
+    localVars.logpost=logpost;
 end
 
 if ~strcmpi(type,'prior'),
     localVars.x=x;
-    localVars.logpost=logpost;
 end
 
 b=0;
@@ -240,6 +261,10 @@ stock_gend=gend;
 stock_data=Y;
 save([DirectoryName '/' M_.fname '_data.mat'],'stock_gend','stock_data');
 
+if strcmpi(type,'gsa')
+    return
+end
+
 if ~isnumeric(options_.parallel),
     leaveSlaveOpen = options_.parallel_info.leaveSlaveOpen;
     if options_.parallel_info.leaveSlaveOpen == 0,
@@ -251,10 +276,10 @@ end
 if options_.smoother
     pm3(endo_nbr,gend,ifil(1),B,'Smoothed variables',...
         '',M_.endo_names(1:M_.orig_endo_nbr, :),'tit_tex',M_.endo_names,...
-        varlist,'SmoothedVariables',[M_.dname '/metropolis'],'_smooth');
+        varlist,'SmoothedVariables',DirectoryName,'_smooth');
     pm3(exo_nbr,gend,ifil(2),B,'Smoothed shocks',...
         '',M_.exo_names,'tit_tex',M_.exo_names,...
-        M_.exo_names,'SmoothedShocks',[M_.dname '/metropolis'],'_inno');
+        M_.exo_names,'SmoothedShocks',DirectoryName,'_inno');
     if nvn
         % needs to  be fixed
         %        pm3(endo_nbr,gend,ifil(3),B,'Smoothed measurement errors',...
@@ -266,20 +291,20 @@ end
 if options_.filtered_vars
     pm3(endo_nbr,gend,ifil(1),B,'Updated Variables',...
         '',varlist,'tit_tex',M_.endo_names,...
-        varlist,'UpdatedVariables',[M_.dname '/metropolis'], ...
+        varlist,'UpdatedVariables',DirectoryName, ...
         '_update');
     pm3(endo_nbr,gend+1,ifil(4),B,'One step ahead forecast',...
         '',varlist,'tit_tex',M_.endo_names,...
-        varlist,'FilteredVariables',[M_.dname '/metropolis'],'_filter_step_ahead');
+        varlist,'FilteredVariables',DirectoryName,'_filter_step_ahead');
 end
 
 if options_.forecast
     pm3(endo_nbr,horizon+maxlag,ifil(6),B,'Forecasted variables (mean)',...
         '',varlist,'tit_tex',M_.endo_names,...
-        varlist,'MeanForecast',[M_.dname '/metropolis'],'_forc_mean');
+        varlist,'MeanForecast',DirectoryName,'_forc_mean');
     pm3(endo_nbr,horizon+maxlag,ifil(6),B,'Forecasted variables (point)',...
         '',varlist,'tit_tex',M_.endo_names,...
-        varlist,'PointForecast',[M_.dname '/metropolis'],'_forc_point');
+        varlist,'PointForecast',DirectoryName,'_forc_point');
 end
 
 
diff --git a/matlab/prior_posterior_statistics_core.m b/matlab/prior_posterior_statistics_core.m
index 7596c82b63676969706087d24e850681b0f564bb..d332da95e115d6b5adfa9fa6a0cd510e93312c78 100644
--- a/matlab/prior_posterior_statistics_core.m
+++ b/matlab/prior_posterior_statistics_core.m
@@ -87,13 +87,15 @@ ifil=myinputs.ifil;
 
 if ~strcmpi(type,'prior'),
     x=myinputs.x;
-    logpost=myinputs.logpost;
+    if strcmpi(type,'posterior'),
+        logpost=myinputs.logpost;
+    end
 end
 if whoiam
     Parallel=myinputs.Parallel;
 end
 
-DirectoryName = CheckPath('metropolis');
+DirectoryName = myinputs.DirectoryName;
 
 RemoteFlag = 0;
 if whoiam
@@ -103,7 +105,7 @@ if whoiam
         end
     end
     ifil=ifil(:,whoiam);
-    waitbarString = ['Please wait... Bayesian (posterior) subdraws (' int2str(fpar) 'of' int2str(B) ')...'];
+    waitbarString = ['Please wait... ',type,' subdraws (' int2str(fpar) 'of' int2str(B) ')...'];
     if Parallel(ThisMatlab).Local,
         waitbarTitle=['Local '];
     else
@@ -114,7 +116,7 @@ else
     if exist('OCTAVE_VERSION')
         diary off;
     else
-        h = waitbar(0,'Taking subdraws...');
+        h = waitbar(0,['Taking ',type,' subdraws...']);
     end
 
 end
@@ -142,7 +144,11 @@ for b=fpar:B
 
     else
         deep = x(b,:);
-        logpo = logpost(b);
+        if strcmpi(type,'posterior')
+            logpo = logpost(b);
+        else
+            logpo = evaluate_posterior_kernel(deep');
+        end
     end
     set_all_parameters(deep);
     [dr,info] = resol(oo_.steady_state,0);