From 2e73856f5a4381c760ab1f3a8d9c423aa0702d46 Mon Sep 17 00:00:00 2001
From: Johannes Pfeifer <jpfeifer@gmx.de>
Date: Thu, 7 Dec 2023 21:15:18 +0100
Subject: [PATCH] GSA and identification: move files to namespace

---
 license.txt                                   | 44 +++++-----
 matlab/{gsa => +gsa}/Morris_Measure_Groups.m  |  0
 matlab/{gsa => +gsa}/Sampling_Function_2.m    |  0
 matlab/{gsa/myboxplot.m => +gsa/boxplot.m}    |  4 +-
 matlab/{gsa => +gsa}/cumplot.m                |  0
 .../log_trans_.m => +gsa/log_transform.m}     |  8 +-
 matlab/{gsa => +gsa}/map_calibration.m        | 16 ++--
 .../map_identification.m}                     | 44 +++++-----
 .../monte_carlo_filtering.m}                  | 50 +++++------
 .../monte_carlo_filtering_analysis.m}         | 18 ++--
 .../monte_carlo_moments.m}                    | 10 +--
 .../prior_draw_gsa.m => +gsa/prior_draw.m}    |  3 +-
 matlab/{gsa => +gsa}/priorcdf.m               |  0
 .../reduced_form_mapping.m}                   | 40 ++++-----
 .../reduced_form_screening.m}                 | 16 ++--
 matlab/{dynare_sensitivity.m => +gsa/run.m}   | 18 ++--
 matlab/{gsa => +gsa}/scatter_analysis.m       |  4 +-
 matlab/{gsa => +gsa}/scatter_mcf.m            |  4 +-
 matlab/{gsa => +gsa}/scatter_plots.m          |  2 +-
 matlab/{gsa => +gsa}/set_shocks_param.m       |  0
 .../{gsa/gsa_skewness.m => +gsa/skewness.m}   |  4 +-
 matlab/{gsa/smirnov.m => +gsa/smirnov_test.m} |  6 +-
 .../stab_map_.m => +gsa/stability_mapping.m}  | 22 ++---
 .../stability_mapping_bivariate.m}            |  4 +-
 .../stability_mapping_univariate.m}           | 16 ++--
 .../stand_.m => +gsa/standardize_columns.m}   |  4 +-
 matlab/{gsa => +gsa}/tcrit.m                  |  0
 matlab/{gsa => +gsa}/teff.m                   |  0
 matlab/{gsa => +gsa}/th_moments.m             |  0
 .../analysis.m}                               | 56 ++++++-------
 .../bruteforce.m}                             |  4 +-
 .../checks.m}                                 | 12 +--
 .../checks_via_subsets.m}                     |  8 +-
 matlab/{ => +identification}/cosn.m           |  2 +-
 .../display.m}                                |  8 +-
 .../get_jacobians.m}                          | 12 +--
 .../numerical_objective.m}                    |  4 +-
 .../plot.m}                                   | 26 +++---
 .../run.m}                                    | 82 +++++++++----------
 .../simulated_moment_uncertainty.m            |  0
 matlab/commutation.m                          |  2 +-
 matlab/duplication.m                          |  2 +-
 matlab/fjaco.m                                |  6 +-
 matlab/get_minimal_state_representation.m     |  2 +-
 matlab/get_perturbation_params_derivs.m       |  2 +-
 matlab/list_of_functions_to_be_cleared.m      |  2 +-
 matlab/pruned_state_space_system.m            |  4 +-
 matlab/set_all_parameters.m                   |  2 +-
 preprocessor                                  |  2 +-
 49 files changed, 288 insertions(+), 287 deletions(-)
 rename matlab/{gsa => +gsa}/Morris_Measure_Groups.m (100%)
 rename matlab/{gsa => +gsa}/Sampling_Function_2.m (100%)
 rename matlab/{gsa/myboxplot.m => +gsa/boxplot.m} (97%)
 rename matlab/{gsa => +gsa}/cumplot.m (100%)
 rename matlab/{gsa/log_trans_.m => +gsa/log_transform.m} (94%)
 rename matlab/{gsa => +gsa}/map_calibration.m (97%)
 rename matlab/{gsa/map_ident_.m => +gsa/map_identification.m} (90%)
 rename matlab/{gsa/filt_mc_.m => +gsa/monte_carlo_filtering.m} (95%)
 rename matlab/{gsa/mcf_analysis.m => +gsa/monte_carlo_filtering_analysis.m} (75%)
 rename matlab/{gsa/mc_moments.m => +gsa/monte_carlo_moments.m} (84%)
 rename matlab/{gsa/prior_draw_gsa.m => +gsa/prior_draw.m} (96%)
 rename matlab/{gsa => +gsa}/priorcdf.m (100%)
 rename matlab/{gsa/redform_map.m => +gsa/reduced_form_mapping.m} (95%)
 rename matlab/{gsa/redform_screen.m => +gsa/reduced_form_screening.m} (92%)
 rename matlab/{dynare_sensitivity.m => +gsa/run.m} (95%)
 rename matlab/{gsa => +gsa}/scatter_analysis.m (89%)
 rename matlab/{gsa => +gsa}/scatter_mcf.m (98%)
 rename matlab/{gsa => +gsa}/scatter_plots.m (99%)
 rename matlab/{gsa => +gsa}/set_shocks_param.m (100%)
 rename matlab/{gsa/gsa_skewness.m => +gsa/skewness.m} (95%)
 rename matlab/{gsa/smirnov.m => +gsa/smirnov_test.m} (94%)
 rename matlab/{gsa/stab_map_.m => +gsa/stability_mapping.m} (95%)
 rename matlab/{gsa/stab_map_2.m => +gsa/stability_mapping_bivariate.m} (96%)
 rename matlab/{gsa/stab_map_1.m => +gsa/stability_mapping_univariate.m} (88%)
 rename matlab/{gsa/stand_.m => +gsa/standardize_columns.m} (93%)
 rename matlab/{gsa => +gsa}/tcrit.m (100%)
 rename matlab/{gsa => +gsa}/teff.m (100%)
 rename matlab/{gsa => +gsa}/th_moments.m (100%)
 rename matlab/{identification_analysis.m => +identification/analysis.m} (94%)
 rename matlab/{ident_bruteforce.m => +identification/bruteforce.m} (98%)
 rename matlab/{identification_checks.m => +identification/checks.m} (95%)
 rename matlab/{identification_checks_via_subsets.m => +identification/checks_via_subsets.m} (98%)
 rename matlab/{ => +identification}/cosn.m (98%)
 rename matlab/{disp_identification.m => +identification/display.m} (98%)
 rename matlab/{get_identification_jacobians.m => +identification/get_jacobians.m} (97%)
 rename matlab/{identification_numerical_objective.m => +identification/numerical_objective.m} (97%)
 rename matlab/{plot_identification.m => +identification/plot.m} (95%)
 rename matlab/{dynare_identification.m => +identification/run.m} (94%)
 rename matlab/{ => +identification}/simulated_moment_uncertainty.m (100%)

diff --git a/license.txt b/license.txt
index d6f1336af2..c1bee92b07 100644
--- a/license.txt
+++ b/license.txt
@@ -161,33 +161,33 @@ Comment: The author gave authorization to redistribute
  Journal of Multivariate Analysis, 2008, vol. 99, issue 3,
  pages 542-554.
 
-Files: matlab/gsa/Morris_Measure_Groups.m
-       matlab/gsa/Sampling_Function_2.m
+Files: matlab/+gsa/Morris_Measure_Groups.m
+       matlab/+gsa/Sampling_Function_2.m
 Copyright: 2005 European Commission
-           2012-2017 Dynare Team
+           2012-2013 Dynare Team
 License: GPL-3+
 Comment: Written by Jessica Cariboni and Francesca Campolongo
  Joint Research Centre, The European Commission,
 
-Files: matlab/gsa/cumplot.m
-       matlab/gsa/filt_mc_.m
-       matlab/gsa/gsa_skewness.m
-       matlab/gsa/log_trans_.m
-       matlab/gsa/map_calibration.m
-       matlab/gsa/map_ident_.m
-       matlab/gsa/mcf_analysis.m
-       matlab/gsa/myboxplot.m
-       matlab/gsa/prior_draw_gsa.m
-       matlab/gsa/redform_map.m
-       matlab/gsa/redform_screen.m
-       matlab/gsa/scatter_mcf.m
-       matlab/gsa/smirnov.m
-       matlab/gsa/stab_map_.m
-       matlab/gsa/stab_map_1.m
-       matlab/gsa/stab_map_2.m
-       matlab/gsa/stand_.m
-       matlab/gsa/tcrit.m
-       matlab/gsa/teff.m
+Files: matlab/+gsa/cumplot.m
+       matlab/+gsa/monte_carlo_filtering.m
+       matlab/+gsa/skewness.m
+       matlab/+gsa/log_trans_.m
+       matlab/+gsa/map_calibration.m
+       matlab/+gsa/map_identification.m
+       matlab/+gsa/monte_carlo_filtering_analysis.m
+       matlab/+gsa/boxplot.m
+       matlab/+gsa/prior_draw.m
+       matlab/+gsa/reduced_form_mapping.m
+       matlab/+gsa/reduced_form_screening.m
+       matlab/+gsa/scatter_mcf.m
+       matlab/+gsa/smirnov_test.m
+       matlab/+gsa/stability_mapping.m
+       matlab/+gsa/stability_mapping_univariate.m
+       matlab/+gsa/stability_mapping_bivariate.m
+       matlab/+gsa/standardize_columns.m
+       matlab/+gsa/tcrit.m
+       matlab/+gsa/teff.m
 Copyright: 2011-2018 European Commission
            2011-2023 Dynare Team
 License: GPL-3+
diff --git a/matlab/gsa/Morris_Measure_Groups.m b/matlab/+gsa/Morris_Measure_Groups.m
similarity index 100%
rename from matlab/gsa/Morris_Measure_Groups.m
rename to matlab/+gsa/Morris_Measure_Groups.m
diff --git a/matlab/gsa/Sampling_Function_2.m b/matlab/+gsa/Sampling_Function_2.m
similarity index 100%
rename from matlab/gsa/Sampling_Function_2.m
rename to matlab/+gsa/Sampling_Function_2.m
diff --git a/matlab/gsa/myboxplot.m b/matlab/+gsa/boxplot.m
similarity index 97%
rename from matlab/gsa/myboxplot.m
rename to matlab/+gsa/boxplot.m
index 4d6cf60d10..f893b7a819 100644
--- a/matlab/gsa/myboxplot.m
+++ b/matlab/+gsa/boxplot.m
@@ -1,5 +1,5 @@
-function sout = myboxplot (data,notched,symbol,vertical,maxwhisker)
-% sout = myboxplot (data,notched,symbol,vertical,maxwhisker)
+function sout = boxplot (data,notched,symbol,vertical,maxwhisker)
+% sout = boxplot (data,notched,symbol,vertical,maxwhisker)
 % Creates a box plot
 
 % Copyright © 2010-2023 Dynare Team
diff --git a/matlab/gsa/cumplot.m b/matlab/+gsa/cumplot.m
similarity index 100%
rename from matlab/gsa/cumplot.m
rename to matlab/+gsa/cumplot.m
diff --git a/matlab/gsa/log_trans_.m b/matlab/+gsa/log_transform.m
similarity index 94%
rename from matlab/gsa/log_trans_.m
rename to matlab/+gsa/log_transform.m
index 3dedb694e8..852ddb1870 100644
--- a/matlab/gsa/log_trans_.m
+++ b/matlab/+gsa/log_transform.m
@@ -1,5 +1,5 @@
-function [yy, xdir, isig, lam]=log_trans_(y0,xdir0,isig,lam)
-% [yy, xdir, isig, lam]=log_trans_(y0,xdir0,isig,lam)
+function [yy, xdir, isig, lam]=log_transform(y0,xdir0,isig,lam)
+% [yy, xdir, isig, lam]=log_transform(y0,xdir0,isig,lam)
 % Conduct automatic log transformation lam(yy/isig+lam)
 % Inputs:
 %   - y0    [double]    series to transform
@@ -56,10 +56,10 @@ end
 if nargin==1
     xdir0='';
 end
-f=@(lam,y)gsa_skewness(log(y+lam));
+f=@(lam,y)gsa.skewness(log(y+lam));
 isig=1;
 if ~(max(y0)<0 || min(y0)>0)
-    if gsa_skewness(y0)<0
+    if gsa.skewness(y0)<0
         isig=-1;
         y0=-y0;
     end
diff --git a/matlab/gsa/map_calibration.m b/matlab/+gsa/map_calibration.m
similarity index 97%
rename from matlab/gsa/map_calibration.m
rename to matlab/+gsa/map_calibration.m
index 44703dc05c..aa68e09008 100644
--- a/matlab/gsa/map_calibration.m
+++ b/matlab/+gsa/map_calibration.m
@@ -229,7 +229,7 @@ if ~isempty(indx_irf)
         if ~options_.nograph && length(time_matrix{plot_indx(ij)})==1
             set(0,'currentfigure',h1),
             subplot(nrow,ncol, plot_indx(ij)),
-            hc = cumplot(mat_irf{ij}(:,ik));
+            hc = gsa.cumplot(mat_irf{ij}(:,ik));
             a=axis;
             delete(hc);
             x1val=max(endo_prior_restrictions.irf{ij,4}(1),a(1));
@@ -237,7 +237,7 @@ if ~isempty(indx_irf)
             hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'b');
             hold all,
             set(hp,'FaceColor', [0.7 0.8 1])
-            hc = cumplot(mat_irf{ij}(:,ik));
+            hc = gsa.cumplot(mat_irf{ij}(:,ik));
             set(hc,'color','k','linewidth',2)
             hold off,
             %         hold off,
@@ -259,7 +259,7 @@ if ~isempty(indx_irf)
         end
         options_mcf.title = atitle0;
         if ~isempty(indx1) && ~isempty(indx2)
-            mcf_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_);
+            gsa.monte_carlo_filtering_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_);
         end
     end
     for ij=1:nbr_irf_couples
@@ -316,7 +316,7 @@ if ~isempty(indx_irf)
 
                 options_mcf.title = atitle0;
                 if ~isempty(indx1) && ~isempty(indx2)
-                    mcf_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_);
+                    gsa.monte_carlo_filtering_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_);
                 end
             end
         end
@@ -434,7 +434,7 @@ if ~isempty(indx_moment)
         if ~options_.nograph && length(time_matrix{plot_indx(ij)})==1
             set(0,'currentfigure',h2);
             subplot(nrow,ncol,plot_indx(ij)),
-            hc = cumplot(mat_moment{ij}(:,ik));
+            hc = gsa.cumplot(mat_moment{ij}(:,ik));
             a=axis; delete(hc),
             %     hist(mat_moment{ij}),
             x1val=max(endo_prior_restrictions.moment{ij,4}(1),a(1));
@@ -442,7 +442,7 @@ if ~isempty(indx_moment)
             hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'b');
             set(hp,'FaceColor', [0.7 0.8 1])
             hold all
-            hc = cumplot(mat_moment{ij}(:,ik));
+            hc = gsa.cumplot(mat_moment{ij}(:,ik));
             set(hc,'color','k','linewidth',2)
             hold off
             title([endo_prior_restrictions.moment{ij,1},' vs ',endo_prior_restrictions.moment{ij,2},'(',leg,')'],'interpreter','none'),
@@ -463,7 +463,7 @@ if ~isempty(indx_moment)
         end
         options_mcf.title = atitle0;
         if ~isempty(indx1) && ~isempty(indx2)
-            mcf_analysis(xmat, indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_);
+            gsa.monte_carlo_filtering_analysis(xmat, indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_);
         end
     end
     for ij=1:nbr_moment_couples
@@ -520,7 +520,7 @@ if ~isempty(indx_moment)
                 end                
                 options_mcf.title = atitle0;
                 if ~isempty(indx1) && ~isempty(indx2)
-                    mcf_analysis(xmat, indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_);
+                    gsa.monte_carlo_filtering_analysis(xmat, indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_);
                 end
             end
         end
diff --git a/matlab/gsa/map_ident_.m b/matlab/+gsa/map_identification.m
similarity index 90%
rename from matlab/gsa/map_ident_.m
rename to matlab/+gsa/map_identification.m
index 2b7194fa2f..3f9f86b9b6 100644
--- a/matlab/gsa/map_ident_.m
+++ b/matlab/+gsa/map_identification.m
@@ -1,5 +1,5 @@
-function map_ident_(OutputDirectoryName,opt_gsa,M_,oo_,options_,estim_params_,bayestopt_)
-% map_ident_(OutputDirectoryName,opt_gsa,M_,oo_,options_,estim_params_,bayestopt_)
+function map_identification(OutputDirectoryName,opt_gsa,M_,oo_,options_,estim_params_,bayestopt_)
+% map_identification(OutputDirectoryName,opt_gsa,M_,oo_,options_,estim_params_,bayestopt_)
 % Inputs
 %  - OutputDirectoryName [string]    name of the output directory
 %  - opt_gsa             [structure]     GSA options structure
@@ -58,16 +58,16 @@ fname_ = M_.fname;
 
 if opt_gsa.load_ident_files==0
     mss = yys(bayestopt_.mfys,:);
-    mss = teff(mss(:,istable),Nsam,istable);
-    yys = teff(yys(dr.order_var,istable),Nsam,istable);
+    mss = gsa.teff(mss(:,istable),Nsam,istable);
+    yys = gsa.teff(yys(dr.order_var,istable),Nsam,istable);
     if exist('T','var')
-        [vdec, cc, ac] = mc_moments(T, lpmatx, dr, M_, options_, estim_params_);
+        [vdec, cc, ac] = gsa.monte_carlo_moments(T, lpmatx, dr, M_, options_, estim_params_);
     else
         return
     end
 
     if opt_gsa.morris==2
-        pdraws = dynare_identification(M_,oo_,options_,bayestopt_,estim_params_,options_.options_ident,[lpmatx lpmat(istable,:)]);
+        pdraws = identification.run(M_,oo_,options_,bayestopt_,estim_params_,options_.options_ident,[lpmatx lpmat(istable,:)]);
         if ~isempty(pdraws) && max(max(abs(pdraws-[lpmatx lpmat(istable,:)])))==0
             disp(['Sample check OK. Largest difference: ', num2str(max(max(abs(pdraws-[lpmatx lpmat(istable,:)]))))]),
             clear pdraws;
@@ -84,7 +84,7 @@ if opt_gsa.load_ident_files==0
             end
             iplo=iplo+1;
             subplot(2,3,iplo)
-            myboxplot(squeeze(vdec(:,j,:))',[],'.',[],10)
+            gsa.boxplot(squeeze(vdec(:,j,:))',[],'.',[],10)
             set(gca,'xticklabel',' ','fontsize',10,'xtick',1:size(options_.varobs,1))
             set(gca,'xlim',[0.5 size(options_.varobs,1)+0.5])
             set(gca,'ylim',[-2 102])
@@ -105,11 +105,11 @@ if opt_gsa.load_ident_files==0
         end
     end
     for j=1:size(cc,1)
-        cc(j,j,:)=stand_(squeeze(log(cc(j,j,:))))./2;
+        cc(j,j,:)=gsa.standardize_columns(squeeze(log(cc(j,j,:))))./2;
     end
-    [vdec, ~, ir_vdec, ic_vdec] = teff(vdec,Nsam,istable);
-    [cc, ~, ir_cc, ic_cc] = teff(cc,Nsam,istable);
-    [ac, ~, ir_ac, ic_ac] = teff(ac,Nsam,istable);
+    [vdec, ~, ir_vdec, ic_vdec] = gsa.teff(vdec,Nsam,istable);
+    [cc, ~, ir_cc, ic_cc] = gsa.teff(cc,Nsam,istable);
+    [ac, ~, ir_ac, ic_ac] = gsa.teff(ac,Nsam,istable);
 
     nc1= size(T,2);
     endo_nbr = M_.endo_nbr;
@@ -123,7 +123,7 @@ if opt_gsa.load_ident_files==0
     [Aa,Bb] = kalman_transition_matrix(dr,iv,ic);
     A = zeros(size(Aa,1),size(Aa,2)+size(Aa,1),length(istable));
     if ~isempty(lpmatx)
-        M_=set_shocks_param(M_,estim_params_,lpmatx(1,:));
+        M_=gsa.set_shocks_param(M_,estim_params_,lpmatx(1,:));
     end
     A(:,:,1)=[Aa, triu(Bb*M_.Sigma_e*Bb')];
     for j=2:length(istable)
@@ -131,14 +131,14 @@ if opt_gsa.load_ident_files==0
         dr.ghu = T(:, (nc1-M_.exo_nbr+1):end, j);
         [Aa,Bb] = kalman_transition_matrix(dr, iv, ic);
         if ~isempty(lpmatx)
-            M_=set_shocks_param(M_,estim_params_,lpmatx(j,:));
+            M_=gsa.set_shocks_param(M_,estim_params_,lpmatx(j,:));
         end
         A(:,:,j)=[Aa, triu(Bb*M_.Sigma_e*Bb')];
     end
     clear T
     clear lpmatx
 
-    [yt, j0]=teff(A,Nsam,istable);
+    [yt, j0]=gsa.teff(A,Nsam,istable);
     yt = [yys yt];
     if opt_gsa.morris==2
         clear TAU A
@@ -155,7 +155,7 @@ if opt_gsa.morris==1
         if opt_gsa.load_ident_files==0
             SAMorris=NaN(npT,3,size(vdec,2));
             for i=1:size(vdec,2)
-                [~, SAMorris(:,:,i)] = Morris_Measure_Groups(npT, [lpmat0 lpmat], vdec(:,i),nliv);
+                [~, SAMorris(:,:,i)] = gsa.Morris_Measure_Groups(npT, [lpmat0 lpmat], vdec(:,i),nliv);
             end
             SAvdec = squeeze(SAMorris(:,1,:))';
             save([OutputDirectoryName,'/',fname_,'_morris_IDE.mat'],'SAvdec','vdec','ir_vdec','ic_vdec')
@@ -164,7 +164,7 @@ if opt_gsa.morris==1
         end
 
         hh_fig = dyn_figure(options_.nodisplay,'name','Screening identification: variance decomposition');
-        myboxplot(SAvdec,[],'.',[],10)
+        gsa.boxplot(SAvdec,[],'.',[],10)
         set(gca,'xticklabel',' ','fontsize',10,'xtick',1:npT)
         set(gca,'xlim',[0.5 npT+0.5])
         ydum = get(gca,'ylim');
@@ -190,7 +190,7 @@ if opt_gsa.morris==1
         ccac = [mss cc ac];
         SAMorris=NaN(npT,3,size(ccac,2));
         for i=1:size(ccac,2)
-            [~, SAMorris(:,:,i)] = Morris_Measure_Groups(npT, [lpmat0 lpmat], [ccac(:,i)],nliv);
+            [~, SAMorris(:,:,i)] = gsa.Morris_Measure_Groups(npT, [lpmat0 lpmat], [ccac(:,i)],nliv);
         end
         SAcc = squeeze(SAMorris(:,1,:))';
         SAcc = SAcc./(max(SAcc,[],2)*ones(1,npT));
@@ -202,7 +202,7 @@ if opt_gsa.morris==1
     end
 
     hh_fig=dyn_figure(options_.nodisplay,'name','Screening identification: theoretical moments');
-    myboxplot(SAcc,[],'.',[],10)
+    gsa.boxplot(SAcc,[],'.',[],10)
     set(gca,'xticklabel',' ','fontsize',10,'xtick',1:npT)
     set(gca,'xlim',[0.5 npT+0.5])
     set(gca,'ylim',[0 1])
@@ -223,7 +223,7 @@ if opt_gsa.morris==1
     if opt_gsa.load_ident_files==0
         SAMorris=NaN(npT,3,j0);
         for j=1:j0
-            [~, SAMorris(:,:,j)] = Morris_Measure_Groups(npT, [lpmat0 lpmat], yt(:,j),nliv);
+            [~, SAMorris(:,:,j)] = gsa.Morris_Measure_Groups(npT, [lpmat0 lpmat], yt(:,j),nliv);
         end
 
         SAM = squeeze(SAMorris(1:end,1,:));
@@ -249,7 +249,7 @@ if opt_gsa.morris==1
         load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAnorm')
     end
     hh_fig=dyn_figure(options_.nodisplay,'name','Screening identification: model');
-    myboxplot(SAnorm',[],'.',[],10)
+    gsa.boxplot(SAnorm',[],'.',[],10)
     set(gca,'xticklabel',' ','fontsize',10,'xtick',1:npT)
     set(gca,'xlim',[0.5 npT+0.5])
     set(gca,'ylim',[0 1])
@@ -297,7 +297,7 @@ else  % main effects analysis
         catch
             EET=[];
         end
-        ccac = stand_([mss cc ac]);
+        ccac = gsa.standardize_columns([mss cc ac]);
         [pcc, dd] = eig(cov(ccac(istable,:)));
         [latent, isort] = sort(-diag(dd));
         latent = -latent;
@@ -314,7 +314,7 @@ else  % main effects analysis
             if itrans==0
                 y0 = ccac(istable,j);
             elseif itrans==1
-                y0 = log_trans_(ccac(istable,j));
+                y0 = gsa.log_transform(ccac(istable,j));
             else
                 y0 = trank(ccac(istable,j));
             end
diff --git a/matlab/gsa/filt_mc_.m b/matlab/+gsa/monte_carlo_filtering.m
similarity index 95%
rename from matlab/gsa/filt_mc_.m
rename to matlab/+gsa/monte_carlo_filtering.m
index b59f6026c3..69b100b4d3 100644
--- a/matlab/gsa/filt_mc_.m
+++ b/matlab/+gsa/monte_carlo_filtering.m
@@ -1,5 +1,5 @@
-function [rmse_MC, ixx] = filt_mc_(OutDir,options_gsa_,dataset_,dataset_info,M_,oo_,options_,bayestopt_,estim_params_)
-% [rmse_MC, ixx] = filt_mc_(OutDir,options_gsa_,dataset_,dataset_info,M_,oo_,options_,bayestopt_,estim_params_
+function [rmse_MC, ixx] = monte_carlo_filtering(OutDir,options_gsa_,dataset_,dataset_info,M_,oo_,options_,bayestopt_,estim_params_)
+% [rmse_MC, ixx] = monte_carlo_filtering(OutDir,options_gsa_,dataset_,dataset_info,M_,oo_,options_,bayestopt_,estim_params_
 % Inputs:
 %  - OutputDirectoryName [string]       name of the output directory
 %  - options_gsa_        [structure]    GSA options
@@ -288,7 +288,7 @@ options_scatter.OutputDirectoryName = OutDir;
 options_scatter.amcf_name = asname;
 options_scatter.amcf_title = atitle;
 options_scatter.title = tmp_title;
-scatter_analysis(r2_MC, x,options_scatter, options_);
+gsa.scatter_analysis(r2_MC, x,options_scatter, options_);
 % end of visual scatter analysis
 
 if ~options_.opt_gsa.ppost && options_.opt_gsa.lik_only
@@ -320,7 +320,7 @@ if ~options_.opt_gsa.ppost && options_.opt_gsa.lik_only
         options_mcf.nobeha_title_latex = 'worse posterior kernel';
     end
 
-    mcf_analysis(x, ipost(1:nfilt), ipost(nfilt+1:end), options_mcf, M_, options_, bayestopt_, estim_params_);
+    gsa.monte_carlo_filtering_analysis(x, ipost(1:nfilt), ipost(nfilt+1:end), options_mcf, M_, options_, bayestopt_, estim_params_);
     if options_.opt_gsa.pprior
         anam = 'rmse_prior_lik';
         atitle = 'RMSE prior: Log Likelihood Kernel';
@@ -338,7 +338,7 @@ if ~options_.opt_gsa.ppost && options_.opt_gsa.lik_only
         options_mcf.nobeha_title_latex = 'worse likelihood';
     end
 
-    mcf_analysis(x, ilik(1:nfilt), ilik(nfilt+1:end), options_mcf, M_, options_, bayestopt_, estim_params_);
+    gsa.monte_carlo_filtering_analysis(x, ilik(1:nfilt), ilik(nfilt+1:end), options_mcf, M_, options_, bayestopt_, estim_params_);
 
 else
     if options_.opt_gsa.ppost
@@ -367,9 +367,9 @@ else
     SS = zeros(npar+nshock, length(vvarvecm));
     for j = 1:npar+nshock
         for i = 1:length(vvarvecm)
-            [~, P] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j), alpha);
-            [H1] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,1);
-            [H2] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,-1);
+            [~, P] = gsa.smirnov_test(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j), alpha);
+            [H1] = gsa.smirnov_test(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,1);
+            [H2] = gsa.smirnov_test(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,-1);
             if H1==0 && H2==0
                 SS(j,i)=1;
             elseif H1==0
@@ -382,7 +382,7 @@ else
         for i = 1:length(vvarvecm)
             for l = 1:length(vvarvecm)
                 if l~=i && PP(j,i)<alpha && PP(j,l)<alpha
-                    [~,P] = smirnov(x(ixx(1:nfilt0(i),i),j),x(ixx(1:nfilt0(l),l),j), alpha);
+                    [~,P] = gsa.smirnov_test(x(ixx(1:nfilt0(i),i),j),x(ixx(1:nfilt0(l),l),j), alpha);
                     PPV(i,l,j) = P;
                 elseif l==i
                     PPV(i,l,j) = PP(j,i);
@@ -407,11 +407,11 @@ else
                 hh_fig=dyn_figure(options_.nodisplay,'name',[temp_name,' ',int2str(ifig)]);
             end
             subplot(3,3,i-9*(ifig-1))
-            h=cumplot(lnprior(ixx(1:nfilt0(i),i)));
+            h=gsa.cumplot(lnprior(ixx(1:nfilt0(i),i)));
             set(h,'color','blue','linewidth',2)
-            hold on, h=cumplot(lnprior);
+            hold on, h=gsa.cumplot(lnprior);
             set(h,'color','k','linewidth',1)
-            h=cumplot(lnprior(ixx(nfilt0(i)+1:end,i)));
+            h=gsa.cumplot(lnprior(ixx(nfilt0(i)+1:end,i)));
             set(h,'color','red','linewidth',2)
             if options_.TeX
                 title(vvarvecm_tex{i},'interpreter','latex')
@@ -459,11 +459,11 @@ else
                 hh_fig = dyn_figure(options_.nodisplay,'Name',[temp_name,' ',int2str(ifig)]);
             end
             subplot(3,3,i-9*(ifig-1))
-            h=cumplot(likelihood(ixx(1:nfilt0(i),i)));
+            h=gsa.cumplot(likelihood(ixx(1:nfilt0(i),i)));
             set(h,'color','blue','linewidth',2)
-            hold on, h=cumplot(likelihood);
+            hold on, h=gsa.cumplot(likelihood);
             set(h,'color','k','linewidth',1)
-            h=cumplot(likelihood(ixx(nfilt0(i)+1:end,i)));
+            h=gsa.cumplot(likelihood(ixx(nfilt0(i)+1:end,i)));
             set(h,'color','red','linewidth',2)
             if options_.TeX
                 title(vvarvecm_tex{i},'interpreter','latex')
@@ -514,11 +514,11 @@ else
                 hh_fig = dyn_figure(options_.nodisplay,'Name',[temp_name,' ',int2str(ifig)]);
             end
             subplot(3,3,i-9*(ifig-1))
-            h=cumplot(logpo2(ixx(1:nfilt0(i),i)));
+            h=gsa.cumplot(logpo2(ixx(1:nfilt0(i),i)));
             set(h,'color','blue','linewidth',2)
-            hold on, h=cumplot(logpo2);
+            hold on, h=gsa.cumplot(logpo2);
             set(h,'color','k','linewidth',1)
-            h=cumplot(logpo2(ixx(nfilt0(i)+1:end,i)));
+            h=gsa.cumplot(logpo2(ixx(nfilt0(i)+1:end,i)));
             set(h,'color','red','linewidth',2)
             if options_.TeX
                 title(vvarvecm_tex{i},'interpreter','latex')
@@ -756,7 +756,7 @@ else
                 options_mcf.nobeha_title_latex = ['worse fit of ' vvarvecm_tex{iy}];
             end
             options_mcf.title = ['the fit of ' vvarvecm{iy}];
-            mcf_analysis(x, ixx(1:nfilt0(iy),iy), ixx(nfilt0(iy)+1:end,iy), options_mcf, M_, options_, bayestopt_, estim_params_);
+            gsa.monte_carlo_filtering_analysis(x, ixx(1:nfilt0(iy),iy), ixx(nfilt0(iy)+1:end,iy), options_mcf, M_, options_, bayestopt_, estim_params_);
         end
         for iy = 1:length(vvarvecm)
             ipar = find(any(squeeze(PPV(iy,:,:))<alpha));
@@ -764,20 +764,20 @@ else
                 hh_fig = dyn_figure(options_.nodisplay,'name',[temp_name,' observed variable ', vvarvecm{iy}]);
                 for j=1+5*(ix-1):min(length(ipar),5*ix)
                     subplot(2,3,j-5*(ix-1))
-                    h0=cumplot(x(:,ipar(j)));
+                    h0=gsa.cumplot(x(:,ipar(j)));
                     set(h0,'color',[0 0 0])
                     hold on,
                     iobs=find(squeeze(PPV(iy,:,ipar(j)))<alpha);
                     for i = 1:length(vvarvecm)
                         if any(iobs==i) || i==iy
-                            h0=cumplot(x(ixx(1:nfilt0(i),i),ipar(j)));
+                            h0=gsa.cumplot(x(ixx(1:nfilt0(i),i),ipar(j)));
                             if ~isoctave
                                 hcmenu = uicontextmenu;
                                 uimenu(hcmenu,'Label',vvarvecm{i});
                                 set(h0,'uicontextmenu',hcmenu)
                             end
                         else
-                            h0=cumplot(x(ixx(1:nfilt0(i),i),ipar(j))*NaN);
+                            h0=gsa.cumplot(x(ixx(1:nfilt0(i),i),ipar(j))*NaN);
                         end
                         set(h0,'color',a00(i,:),'linewidth',2)
                     end
@@ -829,15 +829,15 @@ else
             hh_fig = dyn_figure(options_.nodisplay,'name',[temp_name,' estimated params and shocks ',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)));
+                h0=gsa.cumplot(x(:,nsnam(j)));
                 set(h0,'color',[0 0 0])
                 hold on,
                 npx=find(SP(nsnam(j),:)==0);
                 for i = 1:length(vvarvecm)
                     if any(npx==i)
-                        h0=cumplot(x(ixx(1:nfilt0(i),i),nsnam(j))*NaN);
+                        h0=gsa.cumplot(x(ixx(1:nfilt0(i),i),nsnam(j))*NaN);
                     else
-                        h0=cumplot(x(ixx(1:nfilt0(i),i),nsnam(j)));
+                        h0=gsa.cumplot(x(ixx(1:nfilt0(i),i),nsnam(j)));
                         if ~isoctave
                             hcmenu = uicontextmenu;
                             uimenu(hcmenu,'Label', vvarvecm{i});
diff --git a/matlab/gsa/mcf_analysis.m b/matlab/+gsa/monte_carlo_filtering_analysis.m
similarity index 75%
rename from matlab/gsa/mcf_analysis.m
rename to matlab/+gsa/monte_carlo_filtering_analysis.m
index 795a664ec8..ed312fcdbe 100644
--- a/matlab/gsa/mcf_analysis.m
+++ b/matlab/+gsa/monte_carlo_filtering_analysis.m
@@ -1,5 +1,5 @@
-function indmcf = mcf_analysis(lpmat, ibeha, inobeha, options_mcf, M_, options_, bayestopt_, estim_params_)
-% indmcf = mcf_analysis(lpmat, ibeha, inobeha, options_mcf, M_, options_, bayestopt_, estim_params_)
+function indmcf = monte_carlo_filtering_analysis(lpmat, ibeha, inobeha, options_mcf, M_, options_, bayestopt_, estim_params_)
+% indmcf = monte_carlo_filtering_analysis(lpmat, ibeha, inobeha, options_mcf, M_, options_, bayestopt_, estim_params_)
 % Inputs:
 % - lpmat               [double]        Monte Carlo matrix
 % - ibeha               [integer]       index of behavioural runs
@@ -66,7 +66,7 @@ if isfield(options_mcf,'xparam1')
 end
 OutputDirectoryName = options_mcf.OutputDirectoryName;
 
-[proba, dproba] = stab_map_1(lpmat, ibeha, inobeha, [],fname_, options_, bayestopt_.name, estim_params_,0);
+[proba, dproba] = gsa.stability_mapping_univariate(lpmat, ibeha, inobeha, [],fname_, options_, bayestopt_.name, estim_params_,0);
 indmcf=find(proba<pvalue_ks);
 [~,jtmp] = sort(proba(indmcf),1,'ascend');
 indmcf = indmcf(jtmp);
@@ -87,11 +87,11 @@ end
 
 if length(ibeha)>10 && length(inobeha)>10
     if options_.TeX
-        indcorr1 = stab_map_2(lpmat(ibeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, beha_title, beha_title_latex);
-        indcorr2 = stab_map_2(lpmat(inobeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, nobeha_title, nobeha_title_latex);
+        indcorr1 = gsa.stability_mapping_bivariate(lpmat(ibeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, beha_title, beha_title_latex);
+        indcorr2 = gsa.stability_mapping_bivariate(lpmat(inobeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, nobeha_title, nobeha_title_latex);
     else
-        indcorr1 = stab_map_2(lpmat(ibeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, beha_title);
-        indcorr2 = stab_map_2(lpmat(inobeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, nobeha_title);
+        indcorr1 = gsa.stability_mapping_bivariate(lpmat(ibeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, beha_title);
+        indcorr2 = gsa.stability_mapping_bivariate(lpmat(inobeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, nobeha_title);
     end    
     indcorr = union(indcorr1(:), indcorr2(:));
     indcorr = indcorr(~ismember(indcorr(:),indmcf));
@@ -104,11 +104,11 @@ if ~isempty(indmcf) && ~options_.nograph
         xx=xparam1(indmcf); 
     end
     if options_.TeX
-        scatter_mcf(lpmat(ibeha,indmcf),lpmat(inobeha,indmcf), param_names_tex(indmcf), ...
+        gsa.scatter_mcf(lpmat(ibeha,indmcf),lpmat(inobeha,indmcf), param_names_tex(indmcf), ...
             '.', [fname_,'_',amcf_name], OutputDirectoryName, amcf_title,xx, options_, ...
             beha_title, nobeha_title, beha_title_latex, nobeha_title_latex)
     else
-        scatter_mcf(lpmat(ibeha,indmcf),lpmat(inobeha,indmcf), param_names_tex(indmcf), ...
+        gsa.scatter_mcf(lpmat(ibeha,indmcf),lpmat(inobeha,indmcf), param_names_tex(indmcf), ...
             '.', [fname_,'_',amcf_name], OutputDirectoryName, amcf_title,xx, options_, ...
             beha_title, nobeha_title)
     end
diff --git a/matlab/gsa/mc_moments.m b/matlab/+gsa/monte_carlo_moments.m
similarity index 84%
rename from matlab/gsa/mc_moments.m
rename to matlab/+gsa/monte_carlo_moments.m
index 31554da12d..510e28e949 100644
--- a/matlab/gsa/mc_moments.m
+++ b/matlab/+gsa/monte_carlo_moments.m
@@ -1,5 +1,5 @@
-function [vdec, cc, ac] = mc_moments(mm, ss, dr, M_, options_, estim_params_)
-% [vdec, cc, ac] = mc_moments(mm, ss, dr, M_, options_,estim_params_)
+function [vdec, cc, ac] = monte_carlo_moments(mm, ss, dr, M_, options_, estim_params_)
+% [vdec, cc, ac] = monte_carlo_moments(mm, ss, dr, M_, options_,estim_params_)
 % Conduct Monte Carlo simulation of second moments for GSA
 % Inputs:
 %  - dr                 [structure]     decision rules
@@ -32,7 +32,7 @@ function [vdec, cc, ac] = mc_moments(mm, ss, dr, M_, options_, estim_params_)
 
 [~, nc1, nsam] = size(mm);
 nobs=length(options_.varobs);
-disp('mc_moments: Computing theoretical moments ...')
+disp('monte_carlo_moments: Computing theoretical moments ...')
 h = dyn_waitbar(0,'Theoretical moments ...');
 vdec = zeros(nobs,M_.exo_nbr,nsam);
 cc = zeros(nobs,nobs,nsam);
@@ -42,9 +42,9 @@ for j=1:nsam
     dr.ghx = mm(:, 1:(nc1-M_.exo_nbr),j);
     dr.ghu = mm(:, (nc1-M_.exo_nbr+1):end, j);
     if ~isempty(ss)
-        M_=set_shocks_param(M_,estim_params_,ss(j,:));
+        M_=gsa.set_shocks_param(M_,estim_params_,ss(j,:));
     end
-    [vdec(:,:,j), corr, autocorr] = th_moments(dr,options_,M_);
+    [vdec(:,:,j), corr, autocorr] = gsa.th_moments(dr,options_,M_);
     cc(:,:,j)=triu(corr);
     dum=NaN(nobs,nobs*options_.ar);
     for i=1:options_.ar
diff --git a/matlab/gsa/prior_draw_gsa.m b/matlab/+gsa/prior_draw.m
similarity index 96%
rename from matlab/gsa/prior_draw_gsa.m
rename to matlab/+gsa/prior_draw.m
index 58731ec0a1..c3b8f8d9df 100644
--- a/matlab/gsa/prior_draw_gsa.m
+++ b/matlab/+gsa/prior_draw.m
@@ -1,4 +1,5 @@
-function pdraw = prior_draw_gsa(M_,bayestopt_,options_,estim_params_,init,rdraw)
+function pdraw = prior_draw(M_,bayestopt_,options_,estim_params_,init,rdraw)
+% pdraw = prior_draw(M_,bayestopt_,options_,estim_params_,init,rdraw)
 % Draws from the prior distributions for use with Sensitivity Toolbox for DYNARE
 %
 % INPUTS
diff --git a/matlab/gsa/priorcdf.m b/matlab/+gsa/priorcdf.m
similarity index 100%
rename from matlab/gsa/priorcdf.m
rename to matlab/+gsa/priorcdf.m
diff --git a/matlab/gsa/redform_map.m b/matlab/+gsa/reduced_form_mapping.m
similarity index 95%
rename from matlab/gsa/redform_map.m
rename to matlab/+gsa/reduced_form_mapping.m
index 8475fd7e12..de575f4cb2 100644
--- a/matlab/gsa/redform_map.m
+++ b/matlab/+gsa/reduced_form_mapping.m
@@ -1,5 +1,5 @@
-function redform_map(dirname,options_gsa_,M_,estim_params_,options_,bayestopt_,oo_)
-% redform_map(dirname,options_gsa_,M_,estim_params_,options_,bayestopt_,oo_)
+function reduced_form_mapping(dirname,options_gsa_,M_,estim_params_,options_,bayestopt_,oo_)
+% reduced_form_mapping(dirname,options_gsa_,M_,estim_params_,options_,bayestopt_,oo_)
 % Inputs:
 % - dirname             [string]    name of the output directory
 % - options_gsa_        [structure] GSA options_
@@ -85,7 +85,7 @@ options_mcf.fname_ = M_.fname;
 options_mcf.OutputDirectoryName = adir;
 
 if ~exist('T','var')
-    stab_map_(dirname,options_gsa_,M_,oo_,options_,bayestopt_,estim_params_);
+    gsa.stability_mapping(dirname,options_gsa_,M_,oo_,options_,bayestopt_,estim_params_);
     if pprior
         load([dirname,filesep,M_.fname,'_prior'],'T');
     else
@@ -182,14 +182,14 @@ for j = 1:length(anamendo)
                         end
                         if ~options_.nograph
                             hf=dyn_figure(options_.nodisplay,'name',['Reduced Form Mapping (Monte Carlo Filtering): ',namendo,' vs ', namexo]);
-                            hc = cumplot(y0);
+                            hc = gsa.cumplot(y0);
                             a=axis; delete(hc);
                             x1val=max(threshold(1),a(1));
                             x2val=min(threshold(2),a(2));
                             hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'b');
                             set(hp,'FaceColor', [0.7 0.8 1])
                             hold all,
-                            hc = cumplot(y0);
+                            hc = gsa.cumplot(y0);
                             set(hc,'color','k','linewidth',2)
                             hold off,
                             if options_.TeX
@@ -218,7 +218,7 @@ for j = 1:length(anamendo)
                         options_mcf.OutputDirectoryName = xdir;
                         if ~isempty(iy) && ~isempty(iyc)
                             fprintf(['%4.1f%% of the ',type,' support matches ',atitle0,'\n'],length(iy)/length(y0)*100)
-                            icheck = mcf_analysis(x0, iy, iyc, options_mcf, M_, options_, bayestopt_, estim_params_);
+                            icheck = gsa.monte_carlo_filtering_analysis(x0, iy, iyc, options_mcf, M_, options_, bayestopt_, estim_params_);
 
                             lpmat=x0(iy,:);
                             if nshocks
@@ -349,14 +349,14 @@ for j = 1:length(anamendo)
                         end
                         if ~options_.nograph
                             hf=dyn_figure(options_.nodisplay,'name',['Reduced Form Mapping (Monte Carlo Filtering): ',namendo,' vs lagged ', namlagendo]);
-                            hc = cumplot(y0);
+                            hc = gsa.cumplot(y0);
                             a=axis; delete(hc);
                             x1val=max(threshold(1),a(1));
                             x2val=min(threshold(2),a(2));
                             hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'b');
                             set(hp,'FaceColor', [0.7 0.8 1])
                             hold all,
-                            hc = cumplot(y0);
+                            hc = gsa.cumplot(y0);
                             set(hc,'color','k','linewidth',2)
                             hold off
                             if options_.TeX
@@ -387,7 +387,7 @@ for j = 1:length(anamendo)
                         if ~isempty(iy) && ~isempty(iyc)
 
                             fprintf(['%4.1f%% of the ',type,' support matches ',atitle0,'\n'],length(iy)/length(y0)*100)
-                            icheck = mcf_analysis(x0, iy, iyc, options_mcf, M_, options_, bayestopt_, estim_params_);
+                            icheck = gsa.monte_carlo_filtering_analysis(x0, iy, iyc, options_mcf, M_, options_, bayestopt_, estim_params_);
 
                             lpmat=x0(iy,:);
                             if nshocks
@@ -476,9 +476,9 @@ end
 if isempty(threshold) && ~options_.nograph
     hh_fig=dyn_figure(options_.nodisplay,'name','Reduced Form GSA');
     if ilog==0
-        myboxplot(si',[],'.',[],10)
+        gsa.boxplot(si',[],'.',[],10)
     else
-        myboxplot(silog',[],'.',[],10)
+        gsa.boxplot(silog',[],'.',[],10)
     end
     xlabel(' ')
     set(gca,'xticklabel',' ','fontsize',10,'xtick',1:np)
@@ -513,7 +513,7 @@ if options_map.prior_range
         x0(:,j)=(x0(:,j)-pd(j,3))./(pd(j,4)-pd(j,3));
     end
 else
-    x0=priorcdf(x0,pshape, pd(:,1), pd(:,2), pd(:,3), pd(:,4));
+    x0=gsa.priorcdf(x0,pshape, pd(:,1), pd(:,2), pd(:,3), pd(:,4));
 end
 
 if ilog
@@ -549,7 +549,7 @@ if iload==0
     ipred = setdiff(1:nrun,ifit);
 
     if ilog
-        [~, ~, isig, lam] = log_trans_(y0(iest));
+        [~, ~, isig, lam] = gsa.log_transform(y0(iest));
         y1 = log(y0*isig+lam);
     end
     if ~options_.nograph
@@ -571,9 +571,9 @@ if iload==0
         title(options_map.title,'interpreter','none')
         subplot(222)
         if ilog
-            hc = cumplot(y1);
+            hc = gsa.cumplot(y1);
         else
-            hc = cumplot(y0);
+            hc = gsa.cumplot(y0);
         end
         set(hc,'color','k','linewidth',2)
         title([options_map.title ' CDF'],'interpreter','none')
@@ -620,7 +620,7 @@ if iload==0
         if nfit<nrun
             if ilog
                 yf = ss_anova_fcast(x0(ipred,:), gsa1);
-                yf = log_trans_(yf,'',isig,lam)+ss_anova_fcast(x0(ipred,:), gsax);
+                yf = gsa.log_transform(yf,'',isig,lam)+ss_anova_fcast(x0(ipred,:), gsax);
             else
                 yf = ss_anova_fcast(x0(ipred,:), gsa_);
             end
@@ -657,7 +657,7 @@ function gsa2 = log2level_map(gsa1, isig, lam)
 nest=length(gsa1.y);
 np = size(gsa1.x0,2);
 gsa2=gsa1;
-gsa2.y = log_trans_(gsa1.y,'',isig,lam);
+gsa2.y = gsa.log_transform(gsa1.y,'',isig,lam);
 gsa2.fit = (exp(gsa1.fit)-lam)*isig;
 gsa2.f0 = mean(gsa2.fit);
 gsa2.out.SSE = sum((gsa2.fit-gsa2.y).^2);
@@ -727,7 +727,7 @@ for jt=1:10
     indy{jt}=find( (y0>post_deciles(jt)) & (y0<=post_deciles(jt+1)));
     leg{jt}=[int2str(jt) '-dec'];
 end
-[proba] = stab_map_1(x0, indy{1}, indy{end}, [], fname, options_, parnames, estim_params_,0);
+[proba] = gsa.stability_mapping_univariate(x0, indy{1}, indy{end}, [], fname, options_, parnames, estim_params_,0);
 indmcf=find(proba<options_mcf.pvalue_ks);
 if isempty(indmcf)
     [~,jtmp] = sort(proba,1,'ascend');
@@ -747,7 +747,7 @@ for jx=1:nbr_par
     subplot(nrow,ncol,jx)
     hold off
     for jt=1:10
-        h=cumplot(x0(indy{jt},indmcf(jx)));
+        h=gsa.cumplot(x0(indy{jt},indmcf(jx)));
         set(h,'color', cmap(jt,:), 'linewidth', 2)
         hold all
     end
@@ -782,7 +782,7 @@ if nargin<5
 end
 if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
     fidTeX = fopen([figpath '.tex'],'w');
-    fprintf(fidTeX,'%% TeX eps-loader file generated by redform_map.m (Dynare).\n');
+    fprintf(fidTeX,'%% TeX eps-loader file generated by reduced_form_mapping.m (Dynare).\n');
     fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
     fprintf(fidTeX,'\\begin{figure}[H]\n');
     fprintf(fidTeX,'\\centering \n');
diff --git a/matlab/gsa/redform_screen.m b/matlab/+gsa/reduced_form_screening.m
similarity index 92%
rename from matlab/gsa/redform_screen.m
rename to matlab/+gsa/reduced_form_screening.m
index 4ef7ad153e..47bf18145f 100644
--- a/matlab/gsa/redform_screen.m
+++ b/matlab/+gsa/reduced_form_screening.m
@@ -1,5 +1,5 @@
-function redform_screen(dirname, options_gsa_, estim_params_, M_, dr, options_, bayestopt_)
-% redform_screen(dirname, options_gsa_, estim_params_, M_, dr, options_, bayestopt_)
+function reduced_form_screening(dirname, options_gsa_, estim_params_, M_, dr, options_, bayestopt_)
+% reduced_form_screening(dirname, options_gsa_, estim_params_, M_, dr, options_, bayestopt_)
 % Conduct reduced form screening
 % Inputs:
 %  - dirname             [string]    name of the output directory
@@ -72,7 +72,7 @@ for j=1:size(anamendo,1)
         namexo_tex = anamexo_tex{jx};
         iexo = strmatch(namexo, M_.exo_names, 'exact');
         if ~isempty(iexo)
-            y0=teff(T(iendo,iexo+nspred,:), kn, istable);
+            y0=gsa.teff(T(iendo,iexo+nspred,:), kn, istable);
             if ~isempty(y0)
                 if mod(iplo,9)==0
                     ifig = ifig+1;
@@ -82,7 +82,7 @@ for j=1:size(anamendo,1)
                 iplo = iplo+1;
                 js = js+1;
                 subplot(3, 3, iplo)
-                [~, SAMorris] = Morris_Measure_Groups(np+nshock, [lpmat0 lpmat], y0, nliv);
+                [~, SAMorris] = gsa.Morris_Measure_Groups(np+nshock, [lpmat0 lpmat], y0, nliv);
                 SAM = squeeze(SAMorris(nshock+1:end,1));
                 SA(:,js) = SAM./(max(SAM)+eps);
                 [~, iso] = sort(-SA(:,js));
@@ -122,7 +122,7 @@ for j=1:size(anamendo,1)
         ilagendo=strmatch(namlagendo, M_.endo_names(dr.order_var(M_.nstatic+1:M_.nstatic+nsok)), 'exact');
 
         if ~isempty(ilagendo)
-            y0=teff(T(iendo,ilagendo,:),kn,istable);
+            y0=gsa.teff(T(iendo,ilagendo,:),kn,istable);
             if ~isempty(y0)
                 if mod(iplo,9)==0
                     ifig=ifig+1;
@@ -132,7 +132,7 @@ for j=1:size(anamendo,1)
                 iplo=iplo+1;
                 js=js+1;
                 subplot(3,3,iplo),
-                [~, SAMorris] = Morris_Measure_Groups(np+nshock, [lpmat0 lpmat], y0,nliv);
+                [~, SAMorris] = gsa.Morris_Measure_Groups(np+nshock, [lpmat0 lpmat], y0,nliv);
                 SAM = squeeze(SAMorris(nshock+1:end,1));
                 SA(:,js)=SAM./(max(SAM)+eps);
                 [~, iso] = sort(-SA(:,js));
@@ -166,7 +166,7 @@ for j=1:size(anamendo,1)
 end
 
 hh_fig=dyn_figure(options_.nodisplay,'Name','Reduced form screening');
-myboxplot(SA',[],'.',[],10)
+gsa.boxplot(SA',[],'.',[],10)
 set(gca,'xticklabel',' ','fontsize',10,'xtick',1:np)
 set(gca,'xlim',[0.5 np+0.5])
 set(gca,'ylim',[0 1])
@@ -191,7 +191,7 @@ if nargin<6
 end
 if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
     fidTeX = fopen([figpath '.tex'],'w');
-    fprintf(fidTeX,'%% TeX eps-loader file generated by redform_screen.m (Dynare).\n');
+    fprintf(fidTeX,'%% TeX eps-loader file generated by reduced_form_screening.m (Dynare).\n');
     fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
     fprintf(fidTeX,'\\begin{figure}[H]\n');
     fprintf(fidTeX,'\\centering \n');
diff --git a/matlab/dynare_sensitivity.m b/matlab/+gsa/run.m
similarity index 95%
rename from matlab/dynare_sensitivity.m
rename to matlab/+gsa/run.m
index 225c19b5c5..c84be3cac0 100644
--- a/matlab/dynare_sensitivity.m
+++ b/matlab/+gsa/run.m
@@ -1,5 +1,5 @@
-function x0=dynare_sensitivity(M_,oo_,options_,bayestopt_,estim_params_,options_gsa)
-% x0=dynare_sensitivity(M_,oo_,options_,bayestopt_,estim_params_,options_gsa)
+function x0=run(M_,oo_,options_,bayestopt_,estim_params_,options_gsa)
+% x0=run(M_,oo_,options_,bayestopt_,estim_params_,options_gsa)
 % Frontend to the Sensitivity Analysis Toolbox for DYNARE
 % Inputs:
 %  - M_                     [structure]     Matlab's structure describing the model
@@ -306,7 +306,7 @@ if (options_gsa.load_stab || options_gsa.load_rmse || options_gsa.load_redform)
 end
 
 if options_gsa.stab && ~options_gsa.ppost
-    x0 = stab_map_(OutputDirectoryName,options_gsa,M_,oo_,options_,bayestopt_,estim_params_);
+    x0 = gsa.stability_mapping(OutputDirectoryName,options_gsa,M_,oo_,options_,bayestopt_,estim_params_);
     if isempty(x0)
         skipline()
         disp('Sensitivity computations stopped: no parameter set provided a unique solution')
@@ -316,11 +316,11 @@ end
 
 options_.opt_gsa = options_gsa;
 if ~isempty(options_gsa.moment_calibration) || ~isempty(options_gsa.irf_calibration)
-    map_calibration(OutputDirectoryName, M_, options_, oo_, estim_params_,bayestopt_);
+    gsa.map_calibration(OutputDirectoryName, M_, options_, oo_, estim_params_,bayestopt_);
 end
 
 if options_gsa.identification
-    map_ident_(OutputDirectoryName,options_gsa,M_,oo_,options_,estim_params_,bayestopt_);
+    gsa.map_identification(OutputDirectoryName,options_gsa,M_,oo_,options_,estim_params_,bayestopt_);
 end
 
 if options_gsa.redform && ~isempty(options_gsa.namendo)
@@ -346,10 +346,10 @@ if options_gsa.redform && ~isempty(options_gsa.namendo)
         save([OutputDirectoryName filesep M_.fname '_mc.mat'],'lpmat','lpmat0','istable','iunstable','iwrong','iindeterm')
         options_gsa.load_stab=1;
 
-        x0 = stab_map_(OutputDirectoryName,options_gsa,M_,oo_,options_,bayestopt_,estim_params_);
+        x0 = gsa.stability_mapping(OutputDirectoryName,options_gsa,M_,oo_,options_,bayestopt_,estim_params_);
     end
     if options_gsa.morris==1
-        redform_screen(OutputDirectoryName,options_gsa, estim_params_, M_, oo_.dr, options_, bayestopt_);
+        gsa.reduced_form_screening(OutputDirectoryName,options_gsa, estim_params_, M_, oo_.dr, options_, bayestopt_);
     else
         % check existence of the SS_ANOVA toolbox
         if isempty(options_gsa.threshold_redform) && ~(exist('gsa_sdp','file')==6 || exist('gsa_sdp','file')==2)
@@ -360,7 +360,7 @@ if options_gsa.redform && ~isempty(options_gsa.namendo)
             fprintf('After obtaining the files, you need to unpack them and set a Matlab Path to those files.\n')
             error('SS-ANOVA-R Toolbox missing!')
         end
-        redform_map(OutputDirectoryName,options_gsa,M_,estim_params_,options_,bayestopt_,oo_);
+        gsa.reduced_form_mapping(OutputDirectoryName,options_gsa,M_,estim_params_,options_,bayestopt_,oo_);
     end
 end
 % RMSE mapping
@@ -415,7 +415,7 @@ if options_gsa.rmse
         end
     end
     clear a;
-    filt_mc_(OutputDirectoryName,options_gsa,dataset_,dataset_info,M_,oo_,options_,bayestopt_,estim_params_);
+    gsa.monte_carlo_filtering(OutputDirectoryName,options_gsa,dataset_,dataset_info,M_,oo_,options_,bayestopt_,estim_params_);
 end
 options_.opt_gsa = options_gsa;
 
diff --git a/matlab/gsa/scatter_analysis.m b/matlab/+gsa/scatter_analysis.m
similarity index 89%
rename from matlab/gsa/scatter_analysis.m
rename to matlab/+gsa/scatter_analysis.m
index 1596271bde..db3cc16267 100644
--- a/matlab/gsa/scatter_analysis.m
+++ b/matlab/+gsa/scatter_analysis.m
@@ -50,8 +50,8 @@ if ~options_.nograph
         xx=xparam1;
     end
     if options_.TeX
-        scatter_plots(lpmat, xdata, param_names_tex, '.', [fname_, '_', amcf_name], OutputDirectoryName, amcf_title, xx, options_)
+        gsa.scatter_plots(lpmat, xdata, param_names_tex, '.', [fname_, '_', amcf_name], OutputDirectoryName, amcf_title, xx, options_)
     else 
-        scatter_plots(lpmat, xdata, param_names, '.', [fname_, '_', amcf_name], OutputDirectoryName, amcf_title, xx, options_)
+        gsa.scatter_plots(lpmat, xdata, param_names, '.', [fname_, '_', amcf_name], OutputDirectoryName, amcf_title, xx, options_)
     end
 end
diff --git a/matlab/gsa/scatter_mcf.m b/matlab/+gsa/scatter_mcf.m
similarity index 98%
rename from matlab/gsa/scatter_mcf.m
rename to matlab/+gsa/scatter_mcf.m
index f10995c798..909734c899 100644
--- a/matlab/gsa/scatter_mcf.m
+++ b/matlab/+gsa/scatter_mcf.m
@@ -96,10 +96,10 @@ for i = 1:p
     for j = 1:p
         h = axes('position',[fL(i),fL(p+1-j),ffl,ffl]);
         if i==j
-            h1=cumplot(X(:,j));
+            h1=gsa.cumplot(X(:,j));
             set(h1,'color',[0 0 1],'LineWidth',1.5)
             hold on,
-            h2=cumplot(Y(:,j));
+            h2=gsa.cumplot(Y(:,j));
             set(h2,'color',[1 0 0],'LineWidth',1.5)
             if ~isempty(xparam1)
                 hold on, plot(xparam1([j j]),[0 1],'k--')
diff --git a/matlab/gsa/scatter_plots.m b/matlab/+gsa/scatter_plots.m
similarity index 99%
rename from matlab/gsa/scatter_plots.m
rename to matlab/+gsa/scatter_plots.m
index 64e76b7150..16e4b126ee 100644
--- a/matlab/gsa/scatter_plots.m
+++ b/matlab/+gsa/scatter_plots.m
@@ -86,7 +86,7 @@ for i = 1:p
     for j = 1:p
         h = axes('position',[fL(i),fL(p+1-j),ffl,ffl]);
         if i==j
-            h1=cumplot(X(:,j));
+            h1=gsa.cumplot(X(:,j));
             set(h,'Tag','cumplot')
             set(h1,'color',[0 0 1],'LineWidth',1.5)
             if ~isempty(xparam1)
diff --git a/matlab/gsa/set_shocks_param.m b/matlab/+gsa/set_shocks_param.m
similarity index 100%
rename from matlab/gsa/set_shocks_param.m
rename to matlab/+gsa/set_shocks_param.m
diff --git a/matlab/gsa/gsa_skewness.m b/matlab/+gsa/skewness.m
similarity index 95%
rename from matlab/gsa/gsa_skewness.m
rename to matlab/+gsa/skewness.m
index 7b6c4d8bf5..a4b75768c2 100644
--- a/matlab/gsa/gsa_skewness.m
+++ b/matlab/+gsa/skewness.m
@@ -1,5 +1,5 @@
-function s=gsa_skewness(y)
-% s=gsa_skewness(y)
+function s=skewness(y)
+% s=skewness(y)
 % Compute normalized skewness of y
 % Inputs:
 %  - y  [double]  input vector
diff --git a/matlab/gsa/smirnov.m b/matlab/+gsa/smirnov_test.m
similarity index 94%
rename from matlab/gsa/smirnov.m
rename to matlab/+gsa/smirnov_test.m
index 0c68141e30..3ca80e3c8b 100644
--- a/matlab/gsa/smirnov.m
+++ b/matlab/+gsa/smirnov_test.m
@@ -1,7 +1,7 @@
-function [H,prob,d] = smirnov(x1 , x2 , alpha, iflag )
+function [H,prob,d] = smirnov_test(x1 , x2 , alpha, iflag )
+% [H,prob,d] = smirnov_test(x1 , x2 , alpha, iflag )
 % Smirnov test for 2 distributions
-%   [H,prob,d] = smirnov(x1 , x2 , alpha, iflag )
-%
+
 % Written by Marco Ratto
 % Joint Research Centre, The European Commission,
 % marco.ratto@ec.europa.eu
diff --git a/matlab/gsa/stab_map_.m b/matlab/+gsa/stability_mapping.m
similarity index 95%
rename from matlab/gsa/stab_map_.m
rename to matlab/+gsa/stability_mapping.m
index 019fd1e4e2..0da8eefb20 100644
--- a/matlab/gsa/stab_map_.m
+++ b/matlab/+gsa/stability_mapping.m
@@ -1,5 +1,5 @@
-function x0 = stab_map_(OutputDirectoryName,opt_gsa,M_,oo_,options_,bayestopt_,estim_params_)
-% x0 = stab_map_(OutputDirectoryName,opt_gsa,M_,oo_,options_,bayestopt_,estim_params_)
+function x0 = stability_mapping(OutputDirectoryName,opt_gsa,M_,oo_,options_,bayestopt_,estim_params_)
+% x0 = stability_mapping(OutputDirectoryName,opt_gsa,M_,oo_,options_,bayestopt_,estim_params_)
 % Mapping of stability regions in the prior ranges applying
 % Monte Carlo filtering techniques.
 %
@@ -37,7 +37,7 @@ function x0 = stab_map_(OutputDirectoryName,opt_gsa,M_,oo_,options_,bayestopt_,e
 % 3) Bivariate plots of significant correlation patterns
 %  ( abs(corrcoef) > alpha2) under the stable and unacceptable subsets
 %
-% USES qmc_sequence, stab_map_1, stab_map_2
+% USES qmc_sequence, gsa.stability_mapping_univariate, gsa.stability_mapping_bivariate
 %
 % Written by Marco Ratto
 % Joint Research Centre, The European Commission,
@@ -147,7 +147,7 @@ if fload==0 %run new MC
     yys=zeros(length(dr_.ys),Nsam);
 
     if opt_gsa.morris == 1
-        [lpmat] = Sampling_Function_2(nliv, np+nshock, ntra, ones(np+nshock, 1), zeros(np+nshock,1), []);
+        [lpmat] = gsa.Sampling_Function_2(nliv, np+nshock, ntra, ones(np+nshock, 1), zeros(np+nshock,1), []);
         lpmat = lpmat.*(nliv-1)/nliv+1/nliv/2;
         Nsam=size(lpmat,1);
         lpmat0 = lpmat(:,1:nshock);
@@ -167,7 +167,7 @@ if fload==0 %run new MC
             end
         end
     end
-    prior_draw_gsa(M_,bayestopt_,options_,estim_params_,1); %initialize
+    gsa.prior_draw(M_,bayestopt_,options_,estim_params_,1); %initialize
     if pprior
         for j=1:nshock
             if opt_gsa.morris~=1
@@ -184,7 +184,7 @@ if fload==0 %run new MC
                 lpmat(:,j)=lpmat(:,j).*(upper_bound-lower_bound)+lower_bound;
             end
         else
-            xx=prior_draw_gsa(M_,bayestopt_,options_,estim_params_,0,[lpmat0 lpmat]);
+            xx=gsa.prior_draw(M_,bayestopt_,options_,estim_params_,0,[lpmat0 lpmat]);
             lpmat0=xx(:,1:nshock);
             lpmat=xx(:,nshock+1:end);
             clear xx;
@@ -500,7 +500,7 @@ if ~isempty(iunstable) || ~isempty(iwrong)
             options_mcf.nobeha_title_latex = 'NO unique Stable Saddle-Path';            
         end
         options_mcf.title = 'unique solution';
-        mcf_analysis(lpmat, istable, itmp, options_mcf, M_, options_, bayestopt_, estim_params_);
+        gsa.monte_carlo_filtering_analysis(lpmat, istable, itmp, options_mcf, M_, options_, bayestopt_, estim_params_);
 
         if ~isempty(iindeterm)
             itmp = isolve(~ismember(isolve,iindeterm));
@@ -513,7 +513,7 @@ if ~isempty(iunstable) || ~isempty(iwrong)
                 options_mcf.nobeha_title_latex = 'indeterminacy';
             end
             options_mcf.title = 'indeterminacy';
-            mcf_analysis(lpmat, itmp, iindeterm, options_mcf, M_, options_, bayestopt_, estim_params_);
+            gsa.monte_carlo_filtering_analysis(lpmat, itmp, iindeterm, options_mcf, M_, options_, bayestopt_, estim_params_);
         end
 
         if ~isempty(ixun)
@@ -527,7 +527,7 @@ if ~isempty(iunstable) || ~isempty(iwrong)
                 options_mcf.nobeha_title_latex = 'explosive solution';
             end
             options_mcf.title = 'instability';
-            mcf_analysis(lpmat, itmp, ixun, options_mcf, M_, options_, bayestopt_, estim_params_);
+            gsa.monte_carlo_filtering_analysis(lpmat, itmp, ixun, options_mcf, M_, options_, bayestopt_, estim_params_);
         end
 
         inorestriction = istable(~ismember(istable,irestriction)); % violation of prior restrictions
@@ -543,7 +543,7 @@ if ~isempty(iunstable) || ~isempty(iwrong)
                 options_mcf.nobeha_title_latex = 'inability to find a solution';
             end
             options_mcf.title = 'inability to find a solution';
-            mcf_analysis(lpmat, itmp, iwrong, options_mcf, M_, options_, bayestopt_, estim_params_);
+            gsa.monte_carlo_filtering_analysis(lpmat, itmp, iwrong, options_mcf, M_, options_, bayestopt_, estim_params_);
         end
 
         if ~isempty(irestriction)
@@ -576,7 +576,7 @@ if ~isempty(iunstable) || ~isempty(iwrong)
                 options_mcf.nobeha_title_latex = 'NO prior IRF/moment calibration';
             end
             options_mcf.title = 'prior restrictions';
-            mcf_analysis([lpmat0 lpmat], irestriction, inorestriction, options_mcf, M_, options_, bayestopt_, estim_params_);
+            gsa.monte_carlo_filtering_analysis([lpmat0 lpmat], irestriction, inorestriction, options_mcf, M_, options_, bayestopt_, estim_params_);
             iok = irestriction(1);
             x0 = [lpmat0(iok,:)'; lpmat(iok,:)'];
         else
diff --git a/matlab/gsa/stab_map_2.m b/matlab/+gsa/stability_mapping_bivariate.m
similarity index 96%
rename from matlab/gsa/stab_map_2.m
rename to matlab/+gsa/stability_mapping_bivariate.m
index c6c06da19f..4d508c6a41 100644
--- a/matlab/gsa/stab_map_2.m
+++ b/matlab/+gsa/stability_mapping_bivariate.m
@@ -1,5 +1,5 @@
-function indcorr = stab_map_2(x,alpha2, pvalue_crit, M_,options_,bayestopt_,estim_params_, case_name_plain, case_name_latex, dirname,xparam1,figtitle,fig_caption_latex)
-% indcorr = stab_map_2(x,alpha2, pvalue_crit, M_,options_,bayestopt_,estim_params_, fnam, fnam_latex, dirname,xparam1,figtitle,fig_caption_latex)
+function indcorr = stability_mapping_bivariate(x,alpha2, pvalue_crit, M_,options_,bayestopt_,estim_params_, case_name_plain, case_name_latex, dirname,xparam1,figtitle,fig_caption_latex)
+% indcorr = stability_mapping_bivariate(x,alpha2, pvalue_crit, M_,options_,bayestopt_,estim_params_, fnam, fnam_latex, dirname,xparam1,figtitle,fig_caption_latex)
 % Inputs:
 %  - x
 %  - alpha2
diff --git a/matlab/gsa/stab_map_1.m b/matlab/+gsa/stability_mapping_univariate.m
similarity index 88%
rename from matlab/gsa/stab_map_1.m
rename to matlab/+gsa/stability_mapping_univariate.m
index d6e6d1680c..56c9d00c95 100644
--- a/matlab/gsa/stab_map_1.m
+++ b/matlab/+gsa/stability_mapping_univariate.m
@@ -1,5 +1,5 @@
-function [proba, dproba] = stab_map_1(lpmat, ibehaviour, inonbehaviour, aname, fname_, options_, parnames, estim_params_, iplot, ipar, dirname, pcrit, atitle)
-% [proba, dproba] = stab_map_1(lpmat, ibehaviour, inonbehaviour, aname, fname_, options_, parnames, estim_params_, iplot, ipar, dirname, pcrit, atitle)
+function [proba, dproba] = stability_mapping_univariate(lpmat, ibehaviour, inonbehaviour, aname, fname_, options_, parnames, estim_params_, iplot, ipar, dirname, pcrit, atitle)
+% [proba, dproba] = stability_mapping_univariate(lpmat, ibehaviour, inonbehaviour, aname, fname_, options_, parnames, estim_params_, iplot, ipar, dirname, pcrit, atitle)
 % Inputs:
 % - lpmat               [double]        Monte Carlo matrix
 % - ibehaviour          [integer]       index of behavioural runs
@@ -18,7 +18,7 @@ function [proba, dproba] = stab_map_1(lpmat, ibehaviour, inonbehaviour, aname, f
 %
 % Plots: dotted lines for BEHAVIOURAL
 %        solid lines for NON BEHAVIOURAL
-% USES smirnov
+% USES gsa.smirnov_test.m
 %
 % Written by Marco Ratto
 % Joint Research Centre, The European Commission,
@@ -71,7 +71,7 @@ end
 proba=NaN(npar,1);
 dproba=NaN(npar,1);
 for j=1:npar
-    [~,P,KSSTAT] = smirnov(lpmat(ibehaviour,j),lpmat(inonbehaviour,j));
+    [~,P,KSSTAT] = gsa.smirnov_test(lpmat(ibehaviour,j),lpmat(inonbehaviour,j));
     proba(j)=P;
     dproba(j)=KSSTAT;
 end
@@ -88,12 +88,12 @@ if iplot && ~options_.nograph
         for j=1+12*(i-1):min(nparplot,12*i)
             subplot(3,4,j-12*(i-1))
             if ~isempty(ibehaviour)
-                h=cumplot(lpmat(ibehaviour,j));
+                h=gsa.cumplot(lpmat(ibehaviour,j));
                 set(h,'color',[0 0 1], 'linestyle',':','LineWidth',1.5)
             end
             hold on
             if ~isempty(inonbehaviour)
-                h=cumplot(lpmat(inonbehaviour,j));
+                h=gsa.cumplot(lpmat(inonbehaviour,j));
                 set(h,'color',[0 0 0],'LineWidth',1.5)
             end
             title([ftit{j},'. p-value ', num2str(proba(ipar(j)),2)],'interpreter','none')
@@ -102,7 +102,7 @@ if iplot && ~options_.nograph
             dyn_saveas(hh_fig,[dirname,filesep,fname_,'_',aname,'_SA_',int2str(i)],options_.nodisplay,options_.graph_format);
             if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
                 fidTeX = fopen([dirname,filesep,fname_,'_',aname,'_SA_',int2str(i) '.tex'],'w');
-                fprintf(fidTeX,'%% TeX eps-loader file generated by stab_map_1.m (Dynare).\n');
+                fprintf(fidTeX,'%% TeX eps-loader file generated by gsa.stability_mapping_univariate.m (Dynare).\n');
                 fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
                 fprintf(fidTeX,'\\begin{figure}[H]\n');
                 fprintf(fidTeX,'\\centering \n');
@@ -117,7 +117,7 @@ if iplot && ~options_.nograph
             dyn_saveas(hh_fig,[dirname,filesep,fname_,'_',aname,'_SA'],options_.nodisplay,options_.graph_format);
             if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
                 fidTeX = fopen([dirname,filesep,fname_,'_',aname,'_SA.tex'],'w');
-                fprintf(fidTeX,'%% TeX eps-loader file generated by stab_map_1.m (Dynare).\n');
+                fprintf(fidTeX,'%% TeX eps-loader file generated by gsa.stability_mapping_univariate.m (Dynare).\n');
                 fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
                 fprintf(fidTeX,'\\begin{figure}[H]\n');
                 fprintf(fidTeX,'\\centering \n');
diff --git a/matlab/gsa/stand_.m b/matlab/+gsa/standardize_columns.m
similarity index 93%
rename from matlab/gsa/stand_.m
rename to matlab/+gsa/standardize_columns.m
index c95ccf1435..9f59f20c68 100644
--- a/matlab/gsa/stand_.m
+++ b/matlab/+gsa/standardize_columns.m
@@ -1,5 +1,5 @@
-function [y, meany, stdy] = stand_(x)
-% [y, meany, stdy] = stand_(x)
+function [y, meany, stdy] = standardize_columns(x)
+% [y, meany, stdy] = standardize_columns(x)
 % Standardise a matrix by columns
 %
 % [x,my,sy]=stand(y)
diff --git a/matlab/gsa/tcrit.m b/matlab/+gsa/tcrit.m
similarity index 100%
rename from matlab/gsa/tcrit.m
rename to matlab/+gsa/tcrit.m
diff --git a/matlab/gsa/teff.m b/matlab/+gsa/teff.m
similarity index 100%
rename from matlab/gsa/teff.m
rename to matlab/+gsa/teff.m
diff --git a/matlab/gsa/th_moments.m b/matlab/+gsa/th_moments.m
similarity index 100%
rename from matlab/gsa/th_moments.m
rename to matlab/+gsa/th_moments.m
diff --git a/matlab/identification_analysis.m b/matlab/+identification/analysis.m
similarity index 94%
rename from matlab/identification_analysis.m
rename to matlab/+identification/analysis.m
index 09c9cefe3a..e2dba28ace 100644
--- a/matlab/identification_analysis.m
+++ b/matlab/+identification/analysis.m
@@ -1,5 +1,5 @@
-function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, derivatives_info, info, error_indicator] = identification_analysis(M_,options_,oo_,bayestopt_,estim_params_,params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, init)
-% [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, derivatives_info, info, error_indicator] = identification_analysis(M_,options_,oo_,bayestopt_,estim_params_,params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, init)
+function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, derivatives_info, info, error_indicator] = analysis(M_,options_,oo_,bayestopt_,estim_params_,params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, init)
+% [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, derivatives_info, info, error_indicator] = analysis(M_,options_,oo_,bayestopt_,estim_params_,params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, init)
 % -------------------------------------------------------------------------
 % This function wraps all identification analysis, i.e. it
 % (1) wraps functions for the theoretical identification analysis based on moments (Iskrev, 2010),
@@ -58,18 +58,18 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide
 %                         indicator on problems
 % -------------------------------------------------------------------------
 % This function is called by
-%   * dynare_identification.m
+%   * identification.run
 % -------------------------------------------------------------------------
 % This function calls
 %   * [M_.fname,'.dynamic']
 %   * dseries
 %   * dsge_likelihood.m
 %   * dyn_vech
-%   * ident_bruteforce
-%   * identification_checks
-%   * identification_checks_via_subsets
+%   * identification.bruteforce
+%   * identification.checks
+%   * identification.checks_via_subsets
 %   * isoctave
-%   * get_identification_jacobians (previously getJJ)
+%   * identification.get_jacobians (previously getJJ)
 %   * matlab_ver_less_than
 %   * prior_bounds
 %   * resol
@@ -120,7 +120,7 @@ if ~isempty(estim_params_)
     M_ = set_all_parameters(params,estim_params_,M_);
 end
 
-%get options (see dynare_identification.m for description of options)
+%get options (see identification.run.m for description of options)
 nlags               = options_ident.ar;
 advanced            = options_ident.advanced;
 replic              = options_ident.replic;
@@ -142,7 +142,7 @@ error_indicator.identification_spectrum=0;
 
 if info(1) == 0 %no errors in solution
     % Compute parameter Jacobians for identification analysis
-    [~, ~, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params_, M_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
+    [~, ~, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = identification.get_jacobians(estim_params_, M_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
     if isempty(dMINIMAL)
         % Komunjer and Ng is not computed if (1) minimality conditions are not fullfilled or (2) there are more shocks and measurement errors than observables, so we need to reset options
         error_indicator.identification_minimal = 1;
@@ -206,7 +206,7 @@ if info(1) == 0 %no errors in solution
                 options_ident_local.no_identification_spectrum = 1; %do not recompute dSPECTRUM
                 options_ident_local.ar = nlags;     %store new lag number
                 options_.ar      = nlags;           %store new lag number
-                [~, ~, ~, ~, ~, ~, MOMENTS, dMOMENTS, ~, ~, ~, ~] = get_identification_jacobians(estim_params_, M_, options_, options_ident_local, indpmodel, indpstderr, indpcorr, indvobs, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
+                [~, ~, ~, ~, ~, ~, MOMENTS, dMOMENTS, ~, ~, ~, ~] = identification.get_jacobians(estim_params_, M_, options_, options_ident_local, indpmodel, indpstderr, indpcorr, indvobs, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
 
                 ind_dMOMENTS = (find(max(abs(dMOMENTS'),[],1) > tol_deriv)); %new index with non-zero rows
             end
@@ -305,7 +305,7 @@ if info(1) == 0 %no errors in solution
                 options_.analytic_derivation = analytic_derivation; %reset option
                 AHess = -AHess; %take negative of hessian
                 if min(eig(AHess))<-tol_rank
-                    error('identification_analysis: Analytic Hessian is not positive semi-definite!')
+                    error('identification.analysis: Analytic Hessian is not positive semi-definite!')
                 end
                 ide_hess.AHess = AHess; %store asymptotic Hessian
                 %normalize asymptotic hessian
@@ -313,9 +313,9 @@ if info(1) == 0 %no errors in solution
                 iflag = any((deltaM.*deltaM)==0); %check if all second-order derivatives wrt to a single parameter are nonzero
                 tildaM = AHess./((deltaM)*(deltaM')); %this normalization is for numerical purposes
                 if iflag || rank(AHess)>rank(tildaM)
-                    [ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(AHess, 0, tol_rank, tol_sv, totparam_nbr);
+                    [ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification.checks(AHess, 0, tol_rank, tol_sv, totparam_nbr);
                 else %use normalized version if possible
-                    [ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(tildaM, 0, tol_rank, tol_sv, totparam_nbr);
+                    [ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification.checks(tildaM, 0, tol_rank, tol_sv, totparam_nbr);
                 end
                 indok = find(max(ide_hess.indno,[],1)==0);
                 ide_uncert_unnormaliz(indok) = sqrt(diag(inv(AHess(indok,indok))))';
@@ -325,7 +325,7 @@ if info(1) == 0 %no errors in solution
                 diag_chh = sum(si_dREDUCEDFORM(:,ind1)'.*temp1)';
                 ind1 = ind1(ind1>stderrparam_nbr+corrparam_nbr);
                 cdynamic = si_dDYNAMIC(:,ind1-stderrparam_nbr-corrparam_nbr)*((AHess(ind1,ind1))\si_dDYNAMIC(:,ind1-stderrparam_nbr-corrparam_nbr)');
-                flag_score = 1; %this is used for the title in plot_identification.m
+                flag_score = 1; %this is used for the title in identification.plot.m
             catch
                 %Asymptotic Hessian via simulation
                 if options_.order > 1       
@@ -336,7 +336,7 @@ if info(1) == 0 %no errors in solution
                     options_.periods                 = periods+100;
                 end                
                 replic = max([replic, length(ind_dMOMENTS)*3]);
-                cmm = simulated_moment_uncertainty(ind_dMOMENTS, periods, replic,options_,M_,oo_); %covariance matrix of moments
+                cmm = identification.simulated_moment_uncertainty(ind_dMOMENTS, periods, replic,options_,M_,oo_); %covariance matrix of moments
                 sd = sqrt(diag(cmm));
                 cc = cmm./(sd*sd');
                 [VV,DD,WW] = eig(cc);
@@ -350,9 +350,9 @@ if info(1) == 0 %no errors in solution
                 iflag = any((deltaM.*deltaM)==0);
                 tildaM = MIM./((deltaM)*(deltaM'));
                 if iflag || rank(MIM)>rank(tildaM)
-                    [ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(MIM, 0, tol_rank, tol_sv, totparam_nbr);
+                    [ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification.checks(MIM, 0, tol_rank, tol_sv, totparam_nbr);
                 else %use normalized version if possible
-                    [ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(tildaM, 0, tol_rank, tol_sv, totparam_nbr);
+                    [ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification.checks(tildaM, 0, tol_rank, tol_sv, totparam_nbr);
                 end
                 indok = find(max(ide_hess.indno,[],1)==0);
                 ind1 = find(ide_hess.ind0);
@@ -363,7 +363,7 @@ if info(1) == 0 %no errors in solution
                 if ~isempty(indok)
                     ide_uncert_unnormaliz(indok) = (sqrt(diag(inv(tildaM(indok,indok))))./deltaM(indok))'; %sqrt(diag(inv(MIM(indok,indok))))';
                 end
-                flag_score = 0; %this is used for the title in plot_identification.m
+                flag_score = 0; %this is used for the title in identification.plot.m
             end % end of computing sample information matrix for identification strength measure
 
             ide_strength_dMOMENTS(indok) = (1./(ide_uncert_unnormaliz(indok)'./abs(params(indok)')));              %this is s_i in Ratto and Iskrev (2011, p.13)
@@ -465,11 +465,11 @@ if info(1) == 0 %no errors in solution
         ide_moments.MOMENTS       = MOMENTS;
 
         if advanced
-            % here we do not normalize (i.e. we set norm_dMOMENTS=1) as the OLS in ident_bruteforce is very sensitive to norm_dMOMENTS
-            [ide_moments.pars, ide_moments.cosndMOMENTS] = ident_bruteforce(M_.dname,M_.fname,dMOMENTS(ind_dMOMENTS,:), max_dim_cova_group, options_.TeX, options_ident.name_tex, options_ident.tittxt, tol_deriv);
+            % here we do not normalize (i.e. we set norm_dMOMENTS=1) as the OLS in identification.bruteforce is very sensitive to norm_dMOMENTS
+            [ide_moments.pars, ide_moments.cosndMOMENTS] = identification.bruteforce(M_.dname,M_.fname,dMOMENTS(ind_dMOMENTS,:), max_dim_cova_group, options_.TeX, options_ident.name_tex, options_ident.tittxt, tol_deriv);
         end
 
-        %here we focus on the unnormalized S and V, which is then used in plot_identification.m and for prior_mc > 1
+        %here we focus on the unnormalized S and V, which is then used in identification.plot.m and for prior_mc > 1
         [~, S, V] = svd(dMOMENTS(ind_dMOMENTS,:),0);
         if size(S,1) == 1
             S = S(1); % edge case that S is not a matrix but a row vector
@@ -522,9 +522,9 @@ if info(1) == 0 %no errors in solution
     
 %% Perform identification checks, i.e. find out which parameters are involved
     if checks_via_subsets
-        % identification_checks_via_subsets is only for debugging
+        % identification.checks_via_subsets is only for debugging
         [ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = ...
-            identification_checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident, error_indicator);
+            identification.checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident, error_indicator);
          if ~error_indicator.identification_minimal
              ide_minimal.minimal_state_space=1;
          else
@@ -532,19 +532,19 @@ if info(1) == 0 %no errors in solution
          end
     else
         [ide_dynamic.cond, ide_dynamic.rank, ide_dynamic.ind0, ide_dynamic.indno, ide_dynamic.ino, ide_dynamic.Mco, ide_dynamic.Pco, ide_dynamic.jweak, ide_dynamic.jweak_pair] = ...
-            identification_checks(dDYNAMIC(ind_dDYNAMIC,:)./norm_dDYNAMIC, 1, tol_rank, tol_sv, modparam_nbr);
+            identification.checks(dDYNAMIC(ind_dDYNAMIC,:)./norm_dDYNAMIC, 1, tol_rank, tol_sv, modparam_nbr);
         if ~options_ident.no_identification_reducedform && ~error_indicator.identification_reducedform
             [ide_reducedform.cond, ide_reducedform.rank, ide_reducedform.ind0, ide_reducedform.indno, ide_reducedform.ino, ide_reducedform.Mco, ide_reducedform.Pco, ide_reducedform.jweak, ide_reducedform.jweak_pair] = ...
-                identification_checks(dREDUCEDFORM(ind_dREDUCEDFORM,:)./norm_dREDUCEDFORM, 1, tol_rank, tol_sv, totparam_nbr);
+                identification.checks(dREDUCEDFORM(ind_dREDUCEDFORM,:)./norm_dREDUCEDFORM, 1, tol_rank, tol_sv, totparam_nbr);
         end
         if ~options_ident.no_identification_moments && ~error_indicator.identification_moments
             [ide_moments.cond, ide_moments.rank, ide_moments.ind0, ide_moments.indno, ide_moments.ino, ide_moments.Mco, ide_moments.Pco, ide_moments.jweak, ide_moments.jweak_pair] = ...
-                identification_checks(dMOMENTS(ind_dMOMENTS,:)./norm_dMOMENTS, 1, tol_rank, tol_sv, totparam_nbr);
+                identification.checks(dMOMENTS(ind_dMOMENTS,:)./norm_dMOMENTS, 1, tol_rank, tol_sv, totparam_nbr);
         end
         if ~options_ident.no_identification_minimal 
             if ~error_indicator.identification_minimal
                 [ide_minimal.cond, ide_minimal.rank, ide_minimal.ind0, ide_minimal.indno, ide_minimal.ino, ide_minimal.Mco, ide_minimal.Pco, ide_minimal.jweak, ide_minimal.jweak_pair] = ...
-                    identification_checks(dMINIMAL(ind_dMINIMAL,:)./norm_dMINIMAL, 2, tol_rank, tol_sv, totparam_nbr);
+                    identification.checks(dMINIMAL(ind_dMINIMAL,:)./norm_dMINIMAL, 2, tol_rank, tol_sv, totparam_nbr);
                 ide_minimal.minimal_state_space=1;
             else
                 ide_minimal.minimal_state_space=0;
@@ -552,7 +552,7 @@ if info(1) == 0 %no errors in solution
         end
         if ~options_ident.no_identification_spectrum && ~error_indicator.identification_spectrum
             [ide_spectrum.cond, ide_spectrum.rank, ide_spectrum.ind0, ide_spectrum.indno, ide_spectrum.ino, ide_spectrum.Mco, ide_spectrum.Pco, ide_spectrum.jweak, ide_spectrum.jweak_pair] = ...
-                identification_checks(tilda_dSPECTRUM, 3, tol_rank, tol_sv, totparam_nbr);
+                identification.checks(tilda_dSPECTRUM, 3, tol_rank, tol_sv, totparam_nbr);
         end
     end
 end
diff --git a/matlab/ident_bruteforce.m b/matlab/+identification/bruteforce.m
similarity index 98%
rename from matlab/ident_bruteforce.m
rename to matlab/+identification/bruteforce.m
index 75229b4e8b..c4be89e401 100644
--- a/matlab/ident_bruteforce.m
+++ b/matlab/+identification/bruteforce.m
@@ -18,7 +18,7 @@ function [pars, cosnJ] = ident_bruteforce(dname,fname,J, max_dim_cova_group, TeX
 %  cosnJ : cosn of each column with the selected group of columns
 % -------------------------------------------------------------------------
 % This function is called by
-%   * identification_analysis.m
+%   * identification.analysis.m
 % =========================================================================
 % Copyright © 2009-2023 Dynare Team
 %
@@ -67,7 +67,7 @@ for ll = 1:max_dim_cova_group
         cosnJ2=zeros(size(tmp2,1),1);
         b=[];
         for jj = 1:size(tmp2,1)
-            [cosnJ2(jj,1), b(:,jj)] = cosn([J(:,ii),J(:,tmp2(jj,:))]);
+            [cosnJ2(jj,1), b(:,jj)] = identification.cosn([J(:,ii),J(:,tmp2(jj,:))]);
         end
         cosnJ(ii,ll) = max(cosnJ2(:,1));
         if cosnJ(ii,ll)>tol_deriv
diff --git a/matlab/identification_checks.m b/matlab/+identification/checks.m
similarity index 95%
rename from matlab/identification_checks.m
rename to matlab/+identification/checks.m
index be54d1be11..71c62c012b 100644
--- a/matlab/identification_checks.m
+++ b/matlab/+identification/checks.m
@@ -1,5 +1,5 @@
-function [condX, rankX, ind0, indno, ixno, Mco, Pco, jweak, jweak_pair] = identification_checks(X, test_flag, tol_rank, tol_sv, param_nbr)
-% function [condX, rankX, ind0, indno, ixno, Mco, Pco, jweak, jweak_pair] = identification_checks(X, test_flag, tol_rank, tol_sv, param_nbr)
+function [condX, rankX, ind0, indno, ixno, Mco, Pco, jweak, jweak_pair] = checks(X, test_flag, tol_rank, tol_sv, param_nbr)
+% function [condX, rankX, ind0, indno, ixno, Mco, Pco, jweak, jweak_pair] = checks(X, test_flag, tol_rank, tol_sv, param_nbr)
 % -------------------------------------------------------------------------
 % Checks rank criteria of identification tests and finds out parameter sets
 % that are not identifiable via the nullspace, pairwise correlation
@@ -24,10 +24,10 @@ function [condX, rankX, ind0, indno, ixno, Mco, Pco, jweak, jweak_pair] = identi
 %    * jweak_pair       [(vech) matrix] gives 1 if a couple parameters has Pco=1 (with tolerance tol_rank)
 % -------------------------------------------------------------------------
 % This function is called by
-%   * identification_analysis.m
+%   * identification.analysis.m
 % -------------------------------------------------------------------------
 % This function calls
-%    * cosn
+%    * identification.cosn
 %    * dyn_vech
 %    * vnorm
 % =========================================================================
@@ -141,7 +141,7 @@ if test_flag == 0 || test_flag == 3 % G is a Gram matrix and hence should be a c
 else
     Mco = NaN(param_nbr,1);
     for ii = 1:size(Xparnonzero,2)
-        Mco(ind1(ii),:) = cosn([Xparnonzero(:,ii) , Xparnonzero(:,find([1:1:size(Xparnonzero,2)]~=ii)), Xrest]);
+        Mco(ind1(ii),:) = identification.cosn([Xparnonzero(:,ii) , Xparnonzero(:,find([1:1:size(Xparnonzero,2)]~=ii)), Xrest]);
     end
 end
 
@@ -176,7 +176,7 @@ if test_flag ~= 0
     for ii = 1:size(Xparnonzero,2)
         Pco(ind1(ii),ind1(ii)) = 1;
         for jj = ii+1:size(Xparnonzero,2)
-            Pco(ind1(ii),ind1(jj)) = cosn([Xparnonzero(:,ii),Xparnonzero(:,jj),Xrest]);
+            Pco(ind1(ii),ind1(jj)) = identification.cosn([Xparnonzero(:,ii),Xparnonzero(:,jj),Xrest]);
             Pco(ind1(jj),ind1(ii)) = Pco(ind1(ii),ind1(jj));
         end
     end
diff --git a/matlab/identification_checks_via_subsets.m b/matlab/+identification/checks_via_subsets.m
similarity index 98%
rename from matlab/identification_checks_via_subsets.m
rename to matlab/+identification/checks_via_subsets.m
index 871b882420..4d75e37367 100644
--- a/matlab/identification_checks_via_subsets.m
+++ b/matlab/+identification/checks_via_subsets.m
@@ -1,5 +1,5 @@
-function [ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = identification_checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident,error_indicator)
-%[ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = identification_checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident,error_indicator)
+function [ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident,error_indicator)
+%[ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident,error_indicator)
 % -------------------------------------------------------------------------
 % Finds problematic sets of paramters via checking the necessary rank condition
 % of the Jacobians for all possible combinations of parameters. The rank is
@@ -50,7 +50,7 @@ function [ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal]
 %   * rank:       [integer] rank of Jacobian
 % -------------------------------------------------------------------------
 % This function is called by
-%   * identification_analysis.m
+%   * identification.analysis.m
 % =========================================================================
 % Copyright © 2019-2021 Dynare Team
 %
@@ -161,7 +161,7 @@ end
 
 % initialize for spectrum criteria
 if ~no_identification_spectrum && ~error_indicator.identification_spectrum
-    dSPECTRUM = ide_spectrum.tilda_dSPECTRUM; %tilda dSPECTRUM is normalized dSPECTRUM matrix in identification_analysis.m
+    dSPECTRUM = ide_spectrum.tilda_dSPECTRUM; %tilda dSPECTRUM is normalized dSPECTRUM matrix in identification.analysis.m
     %alternative normalization
     %dSPECTRUM = ide_spectrum.dSPECTRUM;
     %dSPECTRUM(ide_spectrum.ind_dSPECTRUM,:) = dSPECTRUM(ide_spectrum.ind_dSPECTRUM,:)./ide_spectrum.norm_dSPECTRUM; %normalize
diff --git a/matlab/cosn.m b/matlab/+identification/cosn.m
similarity index 98%
rename from matlab/cosn.m
rename to matlab/+identification/cosn.m
index 7ccd1b5bec..a662c245e7 100644
--- a/matlab/cosn.m
+++ b/matlab/+identification/cosn.m
@@ -17,7 +17,7 @@ function [co, b, yhat] = cosn(H)
 %   * y     [n by 1] predicted endogenous values given ols estimation
 % -------------------------------------------------------------------------
 % This function is called by
-%   * identification_checks.m
+%   * identification.checks.m
 %   * ident_bruteforce.m
 % =========================================================================
 % Copyright © 2008-2019 Dynare Team
diff --git a/matlab/disp_identification.m b/matlab/+identification/display.m
similarity index 98%
rename from matlab/disp_identification.m
rename to matlab/+identification/display.m
index c723874329..a0726b8682 100644
--- a/matlab/disp_identification.m
+++ b/matlab/+identification/display.m
@@ -1,5 +1,5 @@
-function disp_identification(pdraws, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, name, options_ident)
-% disp_identification(pdraws, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, name, options_ident)
+function display(pdraws, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, name, options_ident)
+% display(pdraws, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, name, options_ident)
 % -------------------------------------------------------------------------
 % This function displays all identification analysis to the command line
 % =========================================================================
@@ -26,7 +26,7 @@ function disp_identification(pdraws, ide_reducedform, ide_moments, ide_spectrum,
 %   * all output is printed on the command line
 % -------------------------------------------------------------------------
 % This function is called by
-%   * dynare_identification.m
+%   * identification.run
 % =========================================================================
 % Copyright © 2010-2021 Dynare Team
 %
@@ -207,7 +207,7 @@ for jide = 1:4
                 end
             end
 
-            %% display problematic parameters computed by identification_checks_via_subsets
+            %% display problematic parameters computed by identification.checks_via_subsets
         elseif checks_via_subsets
             if ide.rank < size(Jacob,2)
                 no_warning_message_display = 0;
diff --git a/matlab/get_identification_jacobians.m b/matlab/+identification/get_jacobians.m
similarity index 97%
rename from matlab/get_identification_jacobians.m
rename to matlab/+identification/get_jacobians.m
index c60c0d1af3..fc1ba1436d 100644
--- a/matlab/get_identification_jacobians.m
+++ b/matlab/+identification/get_jacobians.m
@@ -1,5 +1,5 @@
-function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
-% [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
+function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_jacobians(estim_params, M_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
+% [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_jacobians(estim_params, M_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
 % previously getJJ.m in Dynare 4.5
 % Sets up the Jacobians needed for identification analysis
 % =========================================================================
@@ -84,7 +84,7 @@ function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dM
 %
 % -------------------------------------------------------------------------
 % This function is called by
-%   * identification_analysis.m
+%   * identification.analysis.m
 % -------------------------------------------------------------------------
 % This function calls
 %   * commutation
@@ -94,7 +94,7 @@ function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dM
 %   * fjaco
 %   * get_perturbation_params_derivs (previously getH)
 %   * get_all_parameters
-%   * identification_numerical_objective (previously thet2tau)
+%   * identification.numerical_objective (previously thet2tau)
 %   * pruned_state_space_system
 %   * vec
 % =========================================================================
@@ -258,7 +258,7 @@ if ~no_identification_moments
     
     if kronflag == -1
         %numerical derivative of autocovariogram
-        dMOMENTS = fjaco(str2func('identification_numerical_objective'), xparam1, 1, estim_params, M_, options_, indpmodel, indpstderr, indvobs, useautocorr, nlags, grid_nbr, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); %[outputflag=1]
+        dMOMENTS = fjaco(str2func('identification.numerical_objective'), xparam1, 1, estim_params, M_, options_, indpmodel, indpstderr, indvobs, useautocorr, nlags, grid_nbr, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); %[outputflag=1]
         dMOMENTS = [dMEAN; dMOMENTS]; %add Jacobian of steady state of VAROBS variables
     else
         dMOMENTS = zeros(obs_nbr + obs_nbr*(obs_nbr+1)/2 + nlags*obs_nbr^2 , totparam_nbr);
@@ -315,7 +315,7 @@ if ~no_identification_spectrum
     IA = eye(size(pruned.A,1));
     if kronflag == -1
         %numerical derivative of spectral density
-        dOmega_tmp = fjaco(str2func('identification_numerical_objective'), xparam1, 2, estim_params, M_, options_, indpmodel, indpstderr, indvobs, useautocorr, nlags, grid_nbr, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); %[outputflag=2]
+        dOmega_tmp = fjaco(str2func('identification.numerical_objective'), xparam1, 2, estim_params, M_, options_, indpmodel, indpstderr, indvobs, useautocorr, nlags, grid_nbr, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); %[outputflag=2]
         kk = 0;
         for ig = 1:length(freqs)
             kk = kk+1;
diff --git a/matlab/identification_numerical_objective.m b/matlab/+identification/numerical_objective.m
similarity index 97%
rename from matlab/identification_numerical_objective.m
rename to matlab/+identification/numerical_objective.m
index 548e687842..ec0ffc7938 100644
--- a/matlab/identification_numerical_objective.m
+++ b/matlab/+identification/numerical_objective.m
@@ -22,7 +22,7 @@ function out = identification_numerical_objective(params, outputflag, estim_para
 % OUTPUTS
 %   out:    dependent on outputflag
 %   *  0:   out = [Yss; vec(A); vec(B); dyn_vech(Sig_e)]; of indvar variables only, in DR order. This is needed to compute dTAU and Komunjer and Ng's D.
-%           Note that Jacobian of Om is computed in get_identification_Jacobians.m (previously getJJ.m) or get_first_order_solution_params_deriv.m (previously getH.m) from Jacobian of B and Sigma_e, because this is more efficient due to some testing with analytical derivatives from An and Schorfheide model
+%           Note that Jacobian of Om is computed in identification.get_jacobians.m (previously getJJ.m) or get_first_order_solution_params_deriv.m (previously getH.m) from Jacobian of B and Sigma_e, because this is more efficient due to some testing with analytical derivatives from An and Schorfheide model
 %   *  1:   out = [vech(cov(Y_t,Y_t)); vec(cov(Y_t,Y_{t-1}); ...; vec(cov(Y_t,Y_{t-nlags})] of indvar variables, in DR order. This is needed to compute Iskrev's J.
 %   *  2:   out = vec(spectral density) with dimension [var_nbr^2*grid_nbr,1] Spectral density of indvar variables evaluated at (grid_nbr/2+1) discretized points in the interval [0;pi]. This is needed for Qu and Tkachenko's G.
 %   * -1:   out = g1(:); of all variables, in DR order. This is needed to compute dLRE.
@@ -32,7 +32,7 @@ function out = identification_numerical_objective(params, outputflag, estim_para
 % Jacobian of the dynamic model equations, and Y_t selected variables
 % -------------------------------------------------------------------------
 % This function is called by
-%   * get_identification_jacobians.m (previously getJJ.m)
+%   * identification.get_jacobians.m (previously getJJ.m)
 % -------------------------------------------------------------------------
 % This function calls
 %   * [M_.fname,'.dynamic']
diff --git a/matlab/plot_identification.m b/matlab/+identification/plot.m
similarity index 95%
rename from matlab/plot_identification.m
rename to matlab/+identification/plot.m
index 035c9a3255..d20531f752 100644
--- a/matlab/plot_identification.m
+++ b/matlab/+identification/plot.m
@@ -1,5 +1,5 @@
-function plot_identification(M_, params, idemoments, idehess, idemodel, idelre, advanced, tittxt, name, IdentifDirectoryName, fname, options_, estim_params_, bayestopt_, tit_TeX, name_tex)
-% plot_identification(M_, params,idemoments,idehess,idemodel, idelre, advanced, tittxt, name, IdentifDirectoryName, fname, options_, estim_params_, bayestopt_, tit_TeX, name_tex)
+function plot(M_, params, idemoments, idehess, idemodel, idelre, advanced, tittxt, name, IdentifDirectoryName, fname, options_, estim_params_, bayestopt_, tit_TeX, name_tex)
+% plot(M_, params,idemoments,idehess,idemodel, idelre, advanced, tittxt, name, IdentifDirectoryName, fname, options_, estim_params_, bayestopt_, tit_TeX, name_tex)
 %
 % INPUTS
 %    o M_                   [structure] model
@@ -156,7 +156,7 @@ if SampleSize == 1
     end
     if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
         fidTeX = fopen([IdentifDirectoryName '/' fname '_ident_strength_' tittxt1,'.tex'],'w');
-        fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n');
+        fprintf(fidTeX,'%% TeX eps-loader file generated by identification.plot.m (Dynare).\n');
         fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
         fprintf(fidTeX,'\\begin{figure}[H]\n');
         fprintf(fidTeX,'\\centering \n');
@@ -203,7 +203,7 @@ if SampleSize == 1
             dyn_saveas(hh_fig,[IdentifDirectoryName '/' fname '_sensitivity_' tittxt1 ],options_.nodisplay,options_.graph_format);
             if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
                 fidTeX = fopen([IdentifDirectoryName '/' fname '_sensitivity_' tittxt1,'.tex'],'w');
-                fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n');
+                fprintf(fidTeX,'%% TeX eps-loader file generated by identification.plot.m (Dynare).\n');
                 fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
                 fprintf(fidTeX,'\\begin{figure}[H]\n');
                 fprintf(fidTeX,'\\centering \n');
@@ -262,7 +262,7 @@ if SampleSize == 1
             dyn_saveas(hh_fig,[ IdentifDirectoryName '/' fname '_ident_collinearity_' tittxt1 '_' int2str(j) ],options_.nodisplay,options_.graph_format);
             if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
                 fidTeX = fopen([ IdentifDirectoryName '/' fname '_ident_collinearity_' tittxt1 '_' int2str(j),'.tex'],'w');
-                fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n');
+                fprintf(fidTeX,'%% TeX eps-loader file generated by identification.plot.m (Dynare).\n');
                 fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
                 fprintf(fidTeX,'\\begin{figure}[H]\n');
                 fprintf(fidTeX,'\\centering \n');
@@ -329,7 +329,7 @@ if SampleSize == 1
         dyn_saveas(f1,[  IdentifDirectoryName '/' fname '_ident_pattern_' tittxt1 '_1' ],options_.nodisplay,options_.graph_format);
         if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
             fidTeX = fopen([  IdentifDirectoryName '/' fname '_ident_pattern_' tittxt1 '_1','.tex'],'w');
-            fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n');
+            fprintf(fidTeX,'%% TeX eps-loader file generated by identification.plot.m (Dynare).\n');
             fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
             fprintf(fidTeX,'\\begin{figure}[H]\n');
             fprintf(fidTeX,'\\centering \n');
@@ -344,7 +344,7 @@ if SampleSize == 1
             dyn_saveas(f2,[  IdentifDirectoryName '/' fname '_ident_pattern_' tittxt1 '_2' ],options_.nodisplay,options_.graph_format);
             if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
                 fidTeX = fopen([  IdentifDirectoryName '/' fname '_ident_pattern_' tittxt1 '_2.tex'],'w');
-                fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n');
+                fprintf(fidTeX,'%% TeX eps-loader file generated by identification.plot.m (Dynare).\n');
                 fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
                 fprintf(fidTeX,'\\begin{figure}[H]\n');
                 fprintf(fidTeX,'\\centering \n');
@@ -392,7 +392,7 @@ else
     dyn_saveas(hh_fig,[ IdentifDirectoryName '/' fname '_MC_sensitivity' ],options_.nodisplay,options_.graph_format);
     if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
         fidTeX = fopen([ IdentifDirectoryName '/' fname '_MC_sensitivity.tex'],'w');
-        fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n');
+        fprintf(fidTeX,'%% TeX eps-loader file generated by identification.plot.m (Dynare).\n');
         fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
         fprintf(fidTeX,'\\begin{figure}[H]\n');
         fprintf(fidTeX,'\\centering \n');
@@ -450,17 +450,17 @@ else
         options_mcf.title = 'MC Highest Condition Number LRE Model';
         ncut=floor(SampleSize/10*9);
         [~,is]=sort(idelre.cond);
-        mcf_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, M_, options_, bayestopt_, estim_params_);
+        gsa.monte_carlo_filtering_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, M_, options_, bayestopt_, estim_params_);
         options_mcf.amcf_name = 'MC_HighestCondNumberModel';
         options_mcf.amcf_title = 'MC Highest Condition Number Model Solution';
         options_mcf.title = 'MC Highest Condition Number Model Solution';
         [~,is]=sort(idemodel.cond);
-        mcf_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, M_, options_, bayestopt_, estim_params_);
+        gsa.monte_carlo_filtering_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, M_, options_, bayestopt_, estim_params_);
         options_mcf.amcf_name = 'MC_HighestCondNumberMoments';
         options_mcf.amcf_title = 'MC Highest Condition Number Model Moments';
         options_mcf.title = 'MC Highest Condition Number Model Moments';
         [~,is]=sort(idemoments.cond);
-        mcf_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, M_, options_, bayestopt_, estim_params_);
+        gsa.monte_carlo_filtering_analysis(params, is(1:ncut), is(ncut+1:end), options_mcf, M_, options_, bayestopt_, estim_params_);
 
         if nparam<5
             f1 = dyn_figure(options_.nodisplay,'Name',[tittxt,' - MC Identification patterns (moments): HIGHEST SV']);
@@ -514,7 +514,7 @@ else
         dyn_saveas(f1,[IdentifDirectoryName '/' fname '_MC_ident_pattern_1' ],options_.nodisplay,options_.graph_format);
         if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
             fidTeX = fopen([IdentifDirectoryName '/' fname '_MC_ident_pattern_1.tex'],'w');
-            fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n');
+            fprintf(fidTeX,'%% TeX eps-loader file generated by identification.plot.m (Dynare).\n');
             fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
             fprintf(fidTeX,'\\begin{figure}[H]\n');
             fprintf(fidTeX,'\\centering \n');
@@ -529,7 +529,7 @@ else
             dyn_saveas(f2,[  IdentifDirectoryName '/' fname '_MC_ident_pattern_2' ],options_.nodisplay,options_.graph_format);
             if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
                 fidTeX = fopen([  IdentifDirectoryName '/' fname '_MC_ident_pattern_2.tex'],'w');
-                fprintf(fidTeX,'%% TeX eps-loader file generated by plot_identification.m (Dynare).\n');
+                fprintf(fidTeX,'%% TeX eps-loader file generated by identification.plot.m (Dynare).\n');
                 fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
                 fprintf(fidTeX,'\\begin{figure}[H]\n');
                 fprintf(fidTeX,'\\centering \n');
diff --git a/matlab/dynare_identification.m b/matlab/+identification/run.m
similarity index 94%
rename from matlab/dynare_identification.m
rename to matlab/+identification/run.m
index 1d9d23dd82..e716fa1ea8 100644
--- a/matlab/dynare_identification.m
+++ b/matlab/+identification/run.m
@@ -1,5 +1,5 @@
-function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, STO_si_dREDUCEDFORM, STO_si_dMOMENTS, STO_dSPECTRUM, STO_dMINIMAL] = dynare_identification(M_,oo_,options_,bayestopt_,estim_params_,options_ident, pdraws0)
-% [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, STO_si_dREDUCEDFORM, STO_si_dMOMENTS, STO_dSPECTRUM, STO_dMINIMAL] = dynare_identification(options_ident, pdraws0)
+function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, STO_si_dREDUCEDFORM, STO_si_dMOMENTS, STO_dSPECTRUM, STO_dMINIMAL] = run(M_,oo_,options_,bayestopt_,estim_params_,options_ident, pdraws0)
+% [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, STO_si_dREDUCEDFORM, STO_si_dMOMENTS, STO_dSPECTRUM, STO_dMINIMAL] = run(options_ident, pdraws0)
 % -------------------------------------------------------------------------
 % This function is called, when the user specifies identification(...); in the mod file. It prepares all identification analysis:
 % (1) set options, local and persistent variables for a new identification
@@ -32,19 +32,19 @@ function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, ST
 % -------------------------------------------------------------------------
 % This function is called by
 %    * driver.m
-%    * map_ident_.m
+%    * gsa.map_identification.m
 % -------------------------------------------------------------------------
 % This function calls
 %    * checkpath
-%    * disp_identification
+%    * identification.display
 %    * dyn_waitbar
 %    * dyn_waitbar_close
 %    * get_all_parameters
 %    * get_posterior_parameters
 %    * get_the_name
-%    * identification_analysis
+%    * identification.analysis
 %    * isoctave
-%    * plot_identification
+%    * identification.plot
 %    * dprior.draw
 %    * set_default_option
 %    * set_prior
@@ -95,7 +95,7 @@ end
 options_ident = set_default_option(options_ident,'gsa_sample_file',0);
     % 0: do not use sample file
     % 1: triggers gsa prior sample
-    % 2: triggers gsa Monte-Carlo sample (i.e. loads a sample corresponding to pprior=0 and ppost=0 in dynare_sensitivity options)
+    % 2: triggers gsa Monte-Carlo sample (i.e. loads a sample corresponding to pprior=0 and ppost=0 in sensitivity.run options)
     % FILENAME: use sample file in provided path
 options_ident = set_default_option(options_ident,'parameter_set','prior_mean');
     % 'calibration':      use values in M_.params and M_.Sigma_e to update estimated stderr, corr and model parameters (get_all_parameters)
@@ -140,7 +140,7 @@ options_ident = set_default_option(options_ident,'tol_rank','robust');
 options_ident = set_default_option(options_ident,'tol_deriv',1.e-8);
     % tolerance level for selecting columns of non-zero derivatives
 options_ident = set_default_option(options_ident,'tol_sv',1.e-3);
-    % tolerance level for selecting non-zero singular values in identification_checks.m
+    % tolerance level for selecting non-zero singular values in identification.checks.m
 options_ident = set_default_option(options_ident,'schur_vec_tol',1e-11);
     % tolerance level used to find nonstationary variables in Schur decomposition of the transition matrix.
 
@@ -181,7 +181,7 @@ if (isfield(options_ident,'no_identification_strength') &&  options_ident.no_ide
     options_ident.no_identification_moments = 0;
 end
 
-%overwrite setting, as dynare_sensitivity does not make use of spectrum and minimal system
+%overwrite setting, as sensitivity.run does not make use of spectrum and minimal system
 if isfield(options_,'opt_gsa') && isfield(options_.opt_gsa,'identification') && options_.opt_gsa.identification == 1
     options_ident.no_identification_minimal = 1;
     options_ident.no_identification_spectrum = 1;
@@ -308,12 +308,12 @@ options_.options_ident = [];
 options_ident = set_default_option(options_ident,'analytic_derivation_mode', options_.analytic_derivation_mode); % if not set by user, inherit default global one
     %  0: efficient sylvester equation method to compute analytical derivatives as in Ratto & Iskrev (2012)
     %  1: kronecker products method to compute analytical derivatives as in Iskrev (2010) (only for order=1)
-    % -1: numerical two-sided finite difference method to compute numerical derivatives of all identification Jacobians using function identification_numerical_objective.m (previously thet2tau.m)
+    % -1: numerical two-sided finite difference method to compute numerical derivatives of all identification Jacobians using function identification.numerical_objective.m (previously thet2tau.m)
     % -2: numerical two-sided finite difference method to compute numerically dYss, dg1, dg2, dg3, d2Yss and d2g1, the identification Jacobians are then computed analytically as with 0
 
 if options_.discretionary_policy || options_.ramsey_policy
     if options_ident.analytic_derivation_mode~=-1
-        fprintf('dynare_identification: discretionary_policy and ramsey_policy require analytic_derivation_mode=-1. Resetting the option.')
+        fprintf('identification.run: discretionary_policy and ramsey_policy require analytic_derivation_mode=-1. Resetting the option.')
         options_ident.analytic_derivation_mode=-1;
     end
 end
@@ -384,7 +384,7 @@ else % no estimated_params block, choose all model parameters and all stderr par
     name_tex = cellfun(@(x) horzcat('$ SE_{', x, '} $'), M_.exo_names_tex, 'UniformOutput', false);
     name_tex = vertcat(name_tex, cellfun(@(x) horzcat('$ ', x, ' $'), M_.param_names_tex, 'UniformOutput', false));
     if ~isequal(M_.H,0)
-        fprintf('\ndynare_identification:: Identification does not support measurement errors (yet) and will ignore them in the following. To test their identifiability, instead define them explicitly as varexo and provide measurement equations in the model definition.\n')
+        fprintf('\nidentification.run:: Identification does not support measurement errors (yet) and will ignore them in the following. To test their identifiability, instead define them explicitly as varexo and provide measurement equations in the model definition.\n')
     end
 end
 options_ident.name_tex = name_tex;
@@ -402,13 +402,13 @@ end
 % settings dependent on number of parameters
 options_ident = set_default_option(options_ident,'max_dim_cova_group',min([2,totparam_nbr-1]));
 options_ident.max_dim_cova_group = min([options_ident.max_dim_cova_group,totparam_nbr-1]);
-    % In brute force search (ident_bruteforce.m) when advanced=1 this option sets the maximum dimension of groups of parameters that best reproduce the behavior of each single model parameter
+    % In brute force search (identification.bruteforce.m) when advanced=1 this option sets the maximum dimension of groups of parameters that best reproduce the behavior of each single model parameter
 
 options_ident = set_default_option(options_ident,'checks_via_subsets',0);
-    % 1: uses identification_checks_via_subsets.m to compute problematic parameter combinations
-    % 0: uses identification_checks.m to compute problematic parameter combinations [default]
+    % 1: uses identification.checks_via_subsets.m to compute problematic parameter combinations
+    % 0: uses identification.checks.m to compute problematic parameter combinations [default]
 options_ident = set_default_option(options_ident,'max_dim_subsets_groups',min([4,totparam_nbr-1]));
-    % In identification_checks_via_subsets.m, when checks_via_subsets=1, this option sets the maximum dimension of groups of parameters for which the corresponding rank criteria is checked
+    % In identification.checks_via_subsets.m, when checks_via_subsets=1, this option sets the maximum dimension of groups of parameters for which the corresponding rank criteria is checked
 
 
 % store identification options
@@ -471,7 +471,7 @@ if iload <=0
     options_ident.tittxt = parameters; %title text for graphs and figures
     % perform identification analysis for single point
     [ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, derivatives_info_point, info, error_indicator_point] = ...
-        identification_analysis(M_,options_,oo_,bayestopt_,estim_params_,params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end implies initialization of persistent variables
+        identification.analysis(M_,options_,oo_,bayestopt_,estim_params_,params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end implies initialization of persistent variables
     if info(1)~=0
         % there are errors in the solution algorithm
         message = get_error_message(info,options_);
@@ -488,7 +488,7 @@ if iload <=0
                 options_ident.tittxt = 'Random_prior_params'; %title text for graphs and figures
                 % perform identification analysis
                 [ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, derivatives_info_point, info, error_indicator_point] = ...
-                    identification_analysis(M_,options_,oo_,bayestopt_,estim_params_,params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1);
+                    identification.analysis(M_,options_,oo_,bayestopt_,estim_params_,params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1);
             end
         end
         if info(1)
@@ -513,10 +513,10 @@ if iload <=0
     save([IdentifDirectoryName '/' fname                '_identif.mat'], 'ide_moments_point', 'ide_spectrum_point', 'ide_minimal_point', 'ide_hess_point', 'ide_reducedform_point', 'ide_dynamic_point', 'store_options_ident');
     save([IdentifDirectoryName '/' fname '_' parameters '_identif.mat'], 'ide_moments_point', 'ide_spectrum_point', 'ide_minimal_point', 'ide_hess_point', 'ide_reducedform_point', 'ide_dynamic_point', 'store_options_ident');
     % display results of identification analysis
-    disp_identification(params, ide_reducedform_point, ide_moments_point, ide_spectrum_point, ide_minimal_point, name, options_ident);
+    identification.display(params, ide_reducedform_point, ide_moments_point, ide_spectrum_point, ide_minimal_point, name, options_ident);
     if ~options_ident.no_identification_strength && ~options_.nograph && ~error_indicator_point.identification_strength && ~error_indicator_point.identification_moments
         % plot (i) identification strength and sensitivity measure based on the moment information matrix and (ii) plot advanced analysis graphs
-        plot_identification(M_,params, ide_moments_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, options_ident.advanced, parameters, name, ...
+        identification.plot(M_,params, ide_moments_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, options_ident.advanced, parameters, name, ...
             IdentifDirectoryName, M_.fname, options_, estim_params_, bayestopt_, parameters_TeX, name_tex);
     end
 
@@ -529,7 +529,7 @@ if iload <=0
         file_index = 0; % initialize counter for files (if options_.MaxNumberOfBytes is reached, we store results in files)
         options_MC = options_ident; %store options structure for Monte Carlo analysis
         options_MC.advanced = 0;    %do not run advanced checking in a Monte Carlo analysis
-        options_ident.checks_via_subsets = 0; % for Monte Carlo analysis currently only identification_checks and not identification_checks_via_subsets is supported
+        options_ident.checks_via_subsets = 0; % for Monte Carlo analysis currently only identification.checks and not identification.checks_via_subsets is supported
     else
         iteration = 1; % iteration equals SampleSize and we are finished
         pdraws = [];   % to have output object otherwise map_ident.m may crash
@@ -543,7 +543,7 @@ if iload <=0
         options_ident.tittxt = []; % clear title text for graphs and figures
         % run identification analysis
         [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, ide_derivatives_info, info, error_indicator] = ...
-            identification_analysis(M_,options_,oo_,bayestopt_,estim_params_,params, indpmodel, indpstderr, indpcorr, options_MC, dataset_info, prior_exist, 0); % the 0 implies that we do not initialize persistent variables anymore
+            identification.analysis(M_,options_,oo_,bayestopt_,estim_params_,params, indpmodel, indpstderr, indpcorr, options_MC, dataset_info, prior_exist, 0); % the 0 implies that we do not initialize persistent variables anymore
 
         if iteration==0 && info(1)==0 % preallocate storage in the first admissable run
             delete([IdentifDirectoryName '/' fname '_identif_*.mat']) % delete previously saved results
@@ -801,25 +801,25 @@ if iload <=0
             end            
             for irun=1:max([maxrun_dDYNAMIC, maxrun_dREDUCEDFORM, maxrun_dMOMENTS, maxrun_dSPECTRUM, maxrun_dMINIMAL])
                 iter=iter+1;
-                % note that this is not the same si_dDYNAMICnorm as computed in identification_analysis
+                % note that this is not the same si_dDYNAMICnorm as computed in identification.analysis
                 % given that we have the MC sample of the Jacobians, we also normalize by the std of the sample of Jacobian entries, to get a fully standardized sensitivity measure
                 si_dDYNAMICnorm(iter,:) = vnorm(STO_si_dDYNAMIC(:,:,irun)./repmat(normalize_STO_DYNAMIC,1,totparam_nbr-(stderrparam_nbr+corrparam_nbr))).*normaliz1((stderrparam_nbr+corrparam_nbr)+1:end);
                 if ~options_MC.no_identification_reducedform && ~isempty(STO_si_dREDUCEDFORM)
-                    % note that this is not the same si_dREDUCEDFORMnorm as computed in identification_analysis
+                    % note that this is not the same si_dREDUCEDFORMnorm as computed in identification.analysis
                     % given that we have the MC sample of the Jacobians, we also normalize by the std of the sample of Jacobian entries, to get a fully standardized sensitivity measure
                     si_dREDUCEDFORMnorm(iter,:) = vnorm(STO_si_dREDUCEDFORM(:,:,irun)./repmat(normalize_STO_REDUCEDFORM,1,totparam_nbr)).*normaliz1;
                 end
                 if ~options_MC.no_identification_moments && ~isempty(STO_si_dMOMENTS)
-                    % note that this is not the same si_dMOMENTSnorm as computed in identification_analysis
+                    % note that this is not the same si_dMOMENTSnorm as computed in identification.analysis
                     % given that we have the MC sample of the Jacobians, we also normalize by the std of the sample of Jacobian entries, to get a fully standardized sensitivity measure
                     si_dMOMENTSnorm(iter,:) = vnorm(STO_si_dMOMENTS(:,:,irun)./repmat(normalize_STO_MOMENTS,1,totparam_nbr)).*normaliz1;
                 end
                 if ~options_MC.no_identification_spectrum && ~isempty(STO_dSPECTRUM)
-                    % note that this is not the same dSPECTRUMnorm as computed in identification_analysis
+                    % note that this is not the same dSPECTRUMnorm as computed in identification.analysis
                     dSPECTRUMnorm(iter,:) = vnorm(STO_dSPECTRUM(:,:,irun)); %not yet used
                 end
                 if ~options_MC.no_identification_minimal && ~isempty(STO_dMINIMAL)
-                    % note that this is not the same dMINIMALnorm as computed in identification_analysis
+                    % note that this is not the same dMINIMALnorm as computed in identification.analysis
                     dMINIMALnorm(iter,:) = vnorm(STO_dMINIMAL(:,:,irun)); %not yet used
                 end
             end
@@ -847,7 +847,7 @@ else
     options_.options_ident      = options_ident;
 end
 
-%% if dynare_identification is called as it own function (not through identification command) and if we load files
+%% if identification.run is called as it own function (not through identification command) and if we load files
 if nargout>3 && iload
     filnam = dir([IdentifDirectoryName '/' fname '_identif_*.mat']);
     STO_si_dDYNAMIC     = [];
@@ -876,10 +876,10 @@ end
 if iload
     %if previous analysis is loaded
     fprintf(['Testing %s\n',parameters]);
-    disp_identification(ide_hess_point.params, ide_reducedform_point, ide_moments_point, ide_spectrum_point, ide_minimal_point, name, options_ident);
+    identification.display(ide_hess_point.params, ide_reducedform_point, ide_moments_point, ide_spectrum_point, ide_minimal_point, name, options_ident);
     if ~options_.nograph && ~error_indicator_point.identification_strength && ~error_indicator_point.identification_moments
         % plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs
-        plot_identification(M_,ide_hess_point.params, ide_moments_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, options_ident.advanced, parameters, name, ...
+        identification.plot(M_,ide_hess_point.params, ide_moments_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, options_ident.advanced, parameters, name, ...
             IdentifDirectoryName, M_.fname, options_, estim_params_, bayestopt_, [], name_tex);
     end
 end
@@ -890,11 +890,11 @@ if SampleSize > 1
     %print results to console but make sure advanced=0
     advanced0 = options_ident.advanced;
     options_ident.advanced = 0;
-    disp_identification(pdraws, IDE_REDUCEDFORM, IDE_MOMENTS, IDE_SPECTRUM, IDE_MINIMAL, name, options_ident);
+    identification.display(pdraws, IDE_REDUCEDFORM, IDE_MOMENTS, IDE_SPECTRUM, IDE_MINIMAL, name, options_ident);
     options_ident.advanced = advanced0; % reset advanced setting
     if ~options_.nograph && isfield(ide_hess_point,'ide_strength_dMOMENTS')
         % plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs
-        plot_identification(M_, pdraws, IDE_MOMENTS, ide_hess_point, IDE_REDUCEDFORM, IDE_DYNAMIC, options_ident.advanced, 'MC sample ', name, ...
+        identification.plot(M_, pdraws, IDE_MOMENTS, ide_hess_point, IDE_REDUCEDFORM, IDE_DYNAMIC, options_ident.advanced, 'MC sample ', name, ...
             IdentifDirectoryName, M_.fname, options_, estim_params_, bayestopt_, [], name_tex);
     end
     %advanced display and plots for MC Sample, i.e. look at draws with highest/lowest condition number
@@ -912,15 +912,15 @@ if SampleSize > 1
                 if ~iload
                     options_ident.tittxt = tittxt; %title text for graphs and figures
                     [ide_moments_max, ide_spectrum_max, ide_minimal_max, ide_hess_max, ide_reducedform_max, ide_dynamic_max, derivatives_info_max, info_max, error_indicator_max] = ...
-                        identification_analysis(M_,options_,oo_,bayestopt_,estim_params_,pdraws(jmax,:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end initializes some persistent variables
+                        identification.analysis(M_,options_,oo_,bayestopt_,estim_params_,pdraws(jmax,:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end initializes some persistent variables
                     save([IdentifDirectoryName '/' fname '_identif.mat'], 'ide_hess_max', 'ide_moments_max', 'ide_spectrum_max', 'ide_minimal_max','ide_reducedform_max', 'ide_dynamic_max', 'jmax', '-append');
                 end
                 advanced0 = options_ident.advanced; options_ident.advanced = 1; % make sure advanced setting is on
-                disp_identification(pdraws(jmax,:), ide_reducedform_max, ide_moments_max, ide_spectrum_max, ide_minimal_max, name, options_ident);
+                identification.display(pdraws(jmax,:), ide_reducedform_max, ide_moments_max, ide_spectrum_max, ide_minimal_max, name, options_ident);
                 options_ident.advanced = advanced0; %reset advanced setting
                 if ~options_.nograph && ~error_indicator_max.identification_strength && ~error_indicator_max.identification_moments
                     % plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs
-                    plot_identification(M_, pdraws(jmax,:), ide_moments_max, ide_hess_max, ide_reducedform_max, ide_dynamic_max, 1, tittxt, name, ...
+                    identification.plot(M_, pdraws(jmax,:), ide_moments_max, ide_hess_max, ide_reducedform_max, ide_dynamic_max, 1, tittxt, name, ...
                         IdentifDirectoryName, M_.fname, options_, estim_params_, bayestopt_, tittxt, name_tex);
                 end
 
@@ -931,15 +931,15 @@ if SampleSize > 1
                 if ~iload
                     options_ident.tittxt = tittxt; %title text for graphs and figures
                     [ide_moments_min, ide_spectrum_min, ide_minimal_min, ide_hess_min, ide_reducedform_min, ide_dynamic_min, ~, ~, error_indicator_min] = ...
-                        identification_analysis(M_,options_,oo_,bayestopt_,estim_params_,pdraws(jmin,:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end initializes persistent variables
+                        identification.analysis(M_,options_,oo_,bayestopt_,estim_params_,pdraws(jmin,:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end initializes persistent variables
                     save([IdentifDirectoryName '/' fname '_identif.mat'], 'ide_hess_min', 'ide_moments_min','ide_spectrum_min','ide_minimal_min','ide_reducedform_min', 'ide_dynamic_min', 'jmin', '-append');
                 end
                 advanced0 = options_ident.advanced; options_ident.advanced = 1; % make sure advanced setting is on
-                disp_identification(pdraws(jmin,:), ide_reducedform_min, ide_moments_min, ide_spectrum_min, ide_minimal_min, name, options_ident);
+                identification.display(pdraws(jmin,:), ide_reducedform_min, ide_moments_min, ide_spectrum_min, ide_minimal_min, name, options_ident);
                 options_ident.advanced = advanced0; %reset advanced setting
                 if ~options_.nograph && ~error_indicator_min.identification_strength && ~error_indicator_min.identification_moments
                     % plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs
-                    plot_identification(M_, pdraws(jmin,:),ide_moments_min,ide_hess_min,ide_reducedform_min,ide_dynamic_min,1,tittxt,name,...
+                    identification.plot(M_, pdraws(jmin,:),ide_moments_min,ide_hess_min,ide_reducedform_min,ide_dynamic_min,1,tittxt,name,...
                         IdentifDirectoryName, M_.fname, options_, estim_params_, bayestopt_, tittxt,name_tex);
                 end
                 % reset nodisplay option
@@ -954,14 +954,14 @@ if SampleSize > 1
                     if ~iload
                         options_ident.tittxt = tittxt; %title text for graphs and figures
                         [ide_moments_(j), ide_spectrum_(j), ide_minimal_(j), ide_hess_(j), ide_reducedform_(j), ide_dynamic_(j), derivatives_info_(j), info_resolve, error_indicator_j] = ...
-                            identification_analysis(M_,options_,oo_,bayestopt_,estim_params_,pdraws(jcrit(j),:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1);
+                            identification.analysis(M_,options_,oo_,bayestopt_,estim_params_,pdraws(jcrit(j),:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1);
                     end
                     advanced0 = options_ident.advanced; options_ident.advanced = 1; %make sure advanced setting is on
-                    disp_identification(pdraws(jcrit(j),:), ide_reducedform_(j), ide_moments_(j), ide_spectrum_(j), ide_minimal_(j), name, options_ident);
+                    identification.display(pdraws(jcrit(j),:), ide_reducedform_(j), ide_moments_(j), ide_spectrum_(j), ide_minimal_(j), name, options_ident);
                     options_ident.advanced = advanced0; % reset advanced
                     if ~options_.nograph && ~error_indicator_j.identification_strength && ~error_indicator_j.identification_moments
                         % plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs
-                        plot_identification(M_, pdraws(jcrit(j),:), ide_moments_(j), ide_hess_(j), ide_reducedform_(j), ide_dynamic_(j), 1, tittxt, name, ...
+                        identification.plot(M_, pdraws(jcrit(j),:), ide_moments_(j), ide_hess_(j), ide_reducedform_(j), ide_dynamic_(j), 1, tittxt, name, ...
                             IdentifDirectoryName, M_.fname, options_, estim_params_, bayestopt_, tittxt, name_tex);
                     end
                 end
diff --git a/matlab/simulated_moment_uncertainty.m b/matlab/+identification/simulated_moment_uncertainty.m
similarity index 100%
rename from matlab/simulated_moment_uncertainty.m
rename to matlab/+identification/simulated_moment_uncertainty.m
diff --git a/matlab/commutation.m b/matlab/commutation.m
index f0c8c6aa5b..1fffe85e8c 100644
--- a/matlab/commutation.m
+++ b/matlab/commutation.m
@@ -14,7 +14,7 @@ function k = commutation(n, m, sparseflag)
 % -------------------------------------------------------------------------
 % This function is called by
 %   * get_perturbation_params_derivs.m (previously getH.m)
-%   * get_identification_jacobians.m (previously getJJ.m)
+%   * identification.get_jacobians.m (previously getJJ.m)
 %   * pruned_state_space_system.m
 % -------------------------------------------------------------------------
 % This function calls
diff --git a/matlab/duplication.m b/matlab/duplication.m
index 9afecf5203..c69a6719f7 100644
--- a/matlab/duplication.m
+++ b/matlab/duplication.m
@@ -11,7 +11,7 @@ function [Dp,DpMPinv] = duplication(p)
 %   DpMPinv: Moore-Penroze inverse of Dp
 % -------------------------------------------------------------------------
 % This function is called by
-%   * get_identification_jacobians.m (previously getJJ.m)
+%   * identification.get_jacobians.m (previously getJJ.m)
 % =========================================================================
 % Copyright © 1997 Tom Minka <minka@microsoft.com>
 % Copyright © 2019 Dynare Team
diff --git a/matlab/fjaco.m b/matlab/fjaco.m
index 3d41787b15..b020ed2698 100644
--- a/matlab/fjaco.m
+++ b/matlab/fjaco.m
@@ -30,7 +30,7 @@ function fjac = fjaco(f,x,varargin)
 ff=feval(f,x,varargin{:});
 
 tol = eps.^(1/3); %some default value
-if strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification_numerical_objective')
+if strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification.numerical_objective')
     tol= varargin{4}.dynatol.x;
 end
 h = tol.*max(abs(x),1);
@@ -40,12 +40,12 @@ fjac = NaN(length(ff),length(x));
 for j=1:length(x)
     xx = x;
     xx(j) = xh1(j); f1=feval(f,xx,varargin{:});
-    if isempty(f1) && (strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification_numerical_objective') )
+    if isempty(f1) && (strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification.numerical_objective') )
         [~,info]=feval(f,xx,varargin{:});
         disp_info_error_identification_perturbation(info,j);
     end
     xx(j) = xh0(j); f0=feval(f,xx,varargin{:});
-    if isempty(f0) && (strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification_numerical_objective') )
+    if isempty(f0) && (strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification.numerical_objective') )
         [~,info]=feval(f,xx,varargin{:});
         disp_info_error_identification_perturbation(info,j)
     end
diff --git a/matlab/get_minimal_state_representation.m b/matlab/get_minimal_state_representation.m
index 277717f3bc..deb6d43bbc 100644
--- a/matlab/get_minimal_state_representation.m
+++ b/matlab/get_minimal_state_representation.m
@@ -53,7 +53,7 @@ function [CheckCO,minns,minSYS] = get_minimal_state_representation(SYS, derivs_f
 %                 Jacobian (wrt to all parameters) of measurement matrix minD
 % -------------------------------------------------------------------------
 % This function is called by
-%   * get_identification_jacobians.m (previously getJJ.m)
+%   * identification.get_jacobians.m (previously getJJ.m)
 % -------------------------------------------------------------------------
 % This function calls
 %   * check_minimality (embedded)
diff --git a/matlab/get_perturbation_params_derivs.m b/matlab/get_perturbation_params_derivs.m
index fab20ef031..6721b3a42a 100644
--- a/matlab/get_perturbation_params_derivs.m
+++ b/matlab/get_perturbation_params_derivs.m
@@ -88,7 +88,7 @@ function DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, dr
 % -------------------------------------------------------------------------
 % This function is called by
 %   * dsge_likelihood.m
-%   * get_identification_jacobians.m
+%   * identification.get_jacobians.m
 % -------------------------------------------------------------------------
 % This function calls
 %   * [fname,'.dynamic']
diff --git a/matlab/list_of_functions_to_be_cleared.m b/matlab/list_of_functions_to_be_cleared.m
index 78a18484ba..200478bb2b 100644
--- a/matlab/list_of_functions_to_be_cleared.m
+++ b/matlab/list_of_functions_to_be_cleared.m
@@ -1,2 +1,2 @@
-list_of_functions = {'discretionary_policy_1', 'dsge_var_likelihood', 'dyn_first_order_solver', 'dyn_waitbar', 'ep_residuals', 'evaluate_likelihood', 'prior_draw_gsa', 'identification_analysis', 'computeDLIK', 'univariate_computeDLIK', 'metropolis_draw', 'flag_implicit_skip_nan', 'mr_hessian', 'masterParallel', 'auxiliary_initialization', 'auxiliary_particle_filter', 'conditional_filter_proposal', 'conditional_particle_filter', 'gaussian_filter', 'gaussian_filter_bank', 'gaussian_mixture_filter', 'gaussian_mixture_filter_bank', 'Kalman_filter', 'online_auxiliary_filter', 'pruned_state_space_system', 'sequential_importance_particle_filter', 'solve_model_for_online_filter', 'prior_draw', 'priordens',...
+list_of_functions = {'discretionary_policy_1', 'dsge_var_likelihood', 'dyn_first_order_solver', 'dyn_waitbar', 'ep_residuals', 'evaluate_likelihood', '+gsa/prior_draw.m', '+identification/analysis.m', 'computeDLIK', 'univariate_computeDLIK', 'metropolis_draw', 'flag_implicit_skip_nan', 'mr_hessian', 'masterParallel', 'auxiliary_initialization', 'auxiliary_particle_filter', 'conditional_filter_proposal', 'conditional_particle_filter', 'gaussian_filter', 'gaussian_filter_bank', 'gaussian_mixture_filter', 'gaussian_mixture_filter_bank', 'Kalman_filter', 'online_auxiliary_filter', 'pruned_state_space_system', 'sequential_importance_particle_filter', 'solve_model_for_online_filter', 'prior_draw', 'priordens',...
     '+occbin/solver.m','+occbin/mkdatap_anticipated_dyn.m','+occbin/mkdatap_anticipated_2constraints_dyn.m','+occbin/match_function.m','+occbin/solve_one_constraint.m','+occbin/solve_two_constraint.m','+occbin/plot/shock_decomposition.m'};
diff --git a/matlab/pruned_state_space_system.m b/matlab/pruned_state_space_system.m
index 3f1f51a143..40d974ca93 100644
--- a/matlab/pruned_state_space_system.m
+++ b/matlab/pruned_state_space_system.m
@@ -80,8 +80,8 @@ function pruned_state_space = pruned_state_space_system(M_, options_, dr, indy,
 %                 parameter Jacobian of E_y
 % -------------------------------------------------------------------------
 % This function is called by
-%   * get_identification_jacobians.m
-%   * identification_numerical_objective.m
+%   * identification.get_jacobians.m
+%   * identification.numerical_objective.m
 % -------------------------------------------------------------------------
 % This function calls
 %   * allVL1.m
diff --git a/matlab/set_all_parameters.m b/matlab/set_all_parameters.m
index c97f98aa26..ed090a1768 100644
--- a/matlab/set_all_parameters.m
+++ b/matlab/set_all_parameters.m
@@ -26,7 +26,7 @@ function M_ = set_all_parameters(xparam1,estim_params_,M_)
 %! @sp 2
 %! @strong{This function is called by:}
 %! @sp 1
-%! @ref{DsgeSmoother}, @ref{dynare_estimation_1}, @ref{@@gsa/filt_mc_}, @ref{identification_analysis}, @ref{PosteriorFilterSmootherAndForecast}, @ref{prior_posterior_statistics_core}, @ref{prior_sampler}
+%! @ref{DsgeSmoother}, @ref{dynare_estimation_1}, @ref{@@gsa.monte_carlo_filtering}, @ref{identification.analysis}, @ref{PosteriorFilterSmootherAndForecast}, @ref{prior_posterior_statistics_core}, @ref{prior_sampler}
 %! @sp 2
 %! @strong{This function calls:}
 %! @sp 2
diff --git a/preprocessor b/preprocessor
index 3dadac8f19..8d0e8cca5c 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit 3dadac8f191dfa1dd660442d5bc4526c4c218149
+Subproject commit 8d0e8cca5cb78b9dde0ecc867ffb0c64d06dd338
-- 
GitLab