diff --git a/matlab/Q6_plication.m b/matlab/+pruned_SS/Q6_plication.m
similarity index 100%
rename from matlab/Q6_plication.m
rename to matlab/+pruned_SS/Q6_plication.m
diff --git a/matlab/allVL1.m b/matlab/+pruned_SS/allVL1.m
similarity index 100%
rename from matlab/allVL1.m
rename to matlab/+pruned_SS/allVL1.m
diff --git a/matlab/bivmom.m b/matlab/+pruned_SS/bivmom.m
similarity index 100%
rename from matlab/bivmom.m
rename to matlab/+pruned_SS/bivmom.m
diff --git a/matlab/prodmom.m b/matlab/+pruned_SS/prodmom.m
similarity index 97%
rename from matlab/prodmom.m
rename to matlab/+pruned_SS/prodmom.m
index f54e541d68dde8b5170046baad526e40986493de..ef42f6fe0b674ddb133cefc64965b38062dcdcc4 100644
--- a/matlab/prodmom.m
+++ b/matlab/+pruned_SS/prodmom.m
@@ -73,7 +73,7 @@ if m==2
         return
     end
    rho = V(1,2)/sqrt(V(1,1)*V(2,2));
-   y = V(1,1)^(nu(1)/2)*V(2,2)^(nu(2)/2)*bivmom(nu,rho);
+   y = V(1,1)^(nu(1)/2)*V(2,2)^(nu(2)/2)*pruned_SS.bivmom(nu,rho);
    return  
 end
 %
diff --git a/matlab/prodmom_deriv.m b/matlab/+pruned_SS/prodmom_deriv.m
similarity index 98%
rename from matlab/prodmom_deriv.m
rename to matlab/+pruned_SS/prodmom_deriv.m
index 2a2271c02b8f8041d9141e63f7a87a04d176f82e..7c7cd69706445790652d02a40b43322f6e48fd87 100644
--- a/matlab/prodmom_deriv.m
+++ b/matlab/+pruned_SS/prodmom_deriv.m
@@ -103,13 +103,13 @@ if m==2
     rho = V(1,2)/sqrt(V(1,1)*V(2,2));
     if nargout > 1
         drho = dC(ii(1),ii(2),:);
-        [tmp,dtmp] = bivmom(nu,rho);
+        [tmp,dtmp] = pruned_SS.bivmom(nu,rho);
         dy = (nu(1)/2)*V(1,1)^(nu(1)/2-1)*dV(1,1,:) * V(2,2)^(nu(2)/2) * tmp...
            + V(1,1)^(nu(1)/2) * (nu(2)/2)*V(2,2)^(nu(2)/2-1)*dV(2,2,:) * tmp...
            + V(1,1)^(nu(1)/2) * V(2,2)^(nu(2)/2) * dtmp * drho;
         dy = reshape(dy,1,size(dV,3));
     else
-        tmp = bivmom(nu,rho);
+        tmp = pruned_SS.bivmom(nu,rho);
     end
     y = V(1,1)^(nu(1)/2)*V(2,2)^(nu(2)/2)*tmp;
     return  
diff --git a/matlab/quadruplication.m b/matlab/+pruned_SS/quadruplication.m
similarity index 100%
rename from matlab/quadruplication.m
rename to matlab/+pruned_SS/quadruplication.m
diff --git a/matlab/pruned_state_space_system.m b/matlab/pruned_state_space_system.m
index bb408bf3f4808d1b83de9d85b4da8a6f1c292fb5..6d71597917c6b11a924752a50858cfb176fab611 100644
--- a/matlab/pruned_state_space_system.m
+++ b/matlab/pruned_state_space_system.m
@@ -424,8 +424,8 @@ if order > 1
     %Compute unique fourth order product moments of u, i.e. unique(E[kron(kron(kron(u,u),u),u)],'stable')
     u_nbr4    = u_nbr*(u_nbr+1)/2*(u_nbr+2)/3*(u_nbr+3)/4;
     if isempty(QPu)
-        QPu       = quadruplication(u_nbr);
-        COMBOS4   = flipud(allVL1(u_nbr, 4)); %all possible (unique) combinations of powers that sum up to four
+        QPu       = pruned_SS.quadruplication(u_nbr);
+        COMBOS4   = flipud(pruned_SS.allVL1(u_nbr, 4)); %all possible (unique) combinations of powers that sum up to four
     end
     E_u_u_u_u = zeros(u_nbr4,1); %only unique entries
     if compute_derivs && (stderrparam_nbr+corrparam_nbr>0)
@@ -433,9 +433,9 @@ if order > 1
     end
     for j4 = 1:size(COMBOS4,1)
         if compute_derivs && (stderrparam_nbr+corrparam_nbr>0)
-            [E_u_u_u_u(j4), dE_u_u_u_u(j4,:)] = prodmom_deriv(E_uu, 1:u_nbr, COMBOS4(j4,:), dE_uu(:,:,1:(stderrparam_nbr+corrparam_nbr)), dr.derivs.dCorrelation_matrix(:,:,1:(stderrparam_nbr+corrparam_nbr)));
+            [E_u_u_u_u(j4), dE_u_u_u_u(j4,:)] = pruned_SS.prodmom_deriv(E_uu, 1:u_nbr, COMBOS4(j4,:), dE_uu(:,:,1:(stderrparam_nbr+corrparam_nbr)), dr.derivs.dCorrelation_matrix(:,:,1:(stderrparam_nbr+corrparam_nbr)));
         else
-            E_u_u_u_u(j4) = prodmom(E_uu, 1:u_nbr, COMBOS4(j4,:));
+            E_u_u_u_u(j4) = pruned_SS.prodmom(E_uu, 1:u_nbr, COMBOS4(j4,:));
         end
     end
     E_xfxf_uu = kron(E_xfxf,E_uu');
@@ -670,8 +670,8 @@ if order > 1
         % Compute unique sixth-order product moments of u, i.e. unique(E[kron(kron(kron(kron(kron(u,u),u),u),u),u)],'stable')
         u_nbr6        = u_nbr*(u_nbr+1)/2*(u_nbr+2)/3*(u_nbr+3)/4*(u_nbr+4)/5*(u_nbr+5)/6;       
         if isempty(Q6Pu)
-            Q6Pu          = Q6_plication(u_nbr);
-            COMBOS6       = flipud(allVL1(u_nbr, 6)); %all possible (unique) combinations of powers that sum up to six
+            Q6Pu          = pruned_SS.Q6_plication(u_nbr);
+            COMBOS6       = flipud(pruned_SS.allVL1(u_nbr, 6)); %all possible (unique) combinations of powers that sum up to six
         end
         E_u_u_u_u_u_u = zeros(u_nbr6,1); %only unique entries
         if compute_derivs && (stderrparam_nbr+corrparam_nbr>0)
@@ -679,9 +679,9 @@ if order > 1
         end
         for j6 = 1:size(COMBOS6,1)
             if compute_derivs && (stderrparam_nbr+corrparam_nbr>0)
-                [E_u_u_u_u_u_u(j6), dE_u_u_u_u_u_u(j6,:)] = prodmom_deriv(E_uu, 1:u_nbr, COMBOS6(j6,:), dE_uu(:,:,1:(stderrparam_nbr+corrparam_nbr)), dr.derivs.dCorrelation_matrix(:,:,1:(stderrparam_nbr+corrparam_nbr)));
+                [E_u_u_u_u_u_u(j6), dE_u_u_u_u_u_u(j6,:)] = pruned_SS.prodmom_deriv(E_uu, 1:u_nbr, COMBOS6(j6,:), dE_uu(:,:,1:(stderrparam_nbr+corrparam_nbr)), dr.derivs.dCorrelation_matrix(:,:,1:(stderrparam_nbr+corrparam_nbr)));
             else
-                E_u_u_u_u_u_u(j6) = prodmom(E_uu, 1:u_nbr, COMBOS6(j6,:));
+                E_u_u_u_u_u_u(j6) = pruned_SS.prodmom(E_uu, 1:u_nbr, COMBOS6(j6,:));
             end
         end
 
diff --git a/matlab/resol.m b/matlab/resol.m
index 0726b309d4cfa874122b40bf8284a1de08002ce8..a99d94a26ac570c54f6e79e193d49b4f56560938 100644
--- a/matlab/resol.m
+++ b/matlab/resol.m
@@ -1,18 +1,18 @@
-function [dr, info, M, oo] = resol(check_flag, M, options, oo)
-
+function [dr, info, M_, oo_] = resol(check_flag, M_, options_, oo_)
+%[dr, info, M_, oo_] = resol(check_flag, M_, options_, oo_)
 % Computes the perturbation based decision rules of the DSGE model (orders 1 to 3)
 %
 % INPUTS
 % - check_flag    [integer]       scalar, equal to 0 if all the approximation is required, equal to 1 if only the eigenvalues are to be computed.
-% - M             [structure]     Matlab's structure describing the model (M_).
-% - options       [structure]     Matlab's structure describing the current options (options_).
-% - oo            [structure]     Matlab's structure containing the results (oo_).
+% - M_            [structure]     Matlab's structure describing the model
+% - options_      [structure]     Matlab's structure describing the current options
+% - oo_           [structure]     Matlab's structure containing the results
 %
 % OUTPUTS
 % - dr            [structure]     Reduced form model.
 % - info          [integer]       scalar or vector, error code.
-% - M             [structure]     Matlab's structure describing the model (M_).
-% - oo            [structure]     Matlab's structure containing the results (oo_).
+% - M_            [structure]     Matlab's structure describing the model
+% - oo_           [structure]     Matlab's structure containing the results
 %
 % REMARKS
 % Possible values for the error codes are:
@@ -49,46 +49,46 @@ function [dr, info, M, oo] = resol(check_flag, M, options, oo)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-if isfield(oo,'dr')
-    if isfield(oo.dr,'kstate')
-        dr.kstate = oo.dr.kstate;
+if isfield(oo_,'dr')
+    if isfield(oo_.dr,'kstate')
+        dr.kstate = oo_.dr.kstate;
     end
-    if isfield(oo.dr,'inv_order_var')
-        dr.inv_order_var = oo.dr.inv_order_var;
+    if isfield(oo_.dr,'inv_order_var')
+        dr.inv_order_var = oo_.dr.inv_order_var;
     end
-    if isfield(oo.dr,'order_var')
-        dr.order_var = oo.dr.order_var;
+    if isfield(oo_.dr,'order_var')
+        dr.order_var = oo_.dr.order_var;
     end
-    if isfield(oo.dr,'restrict_var_list')
-        dr.restrict_var_list = oo.dr.restrict_var_list;
+    if isfield(oo_.dr,'restrict_var_list')
+        dr.restrict_var_list = oo_.dr.restrict_var_list;
     end
-    if isfield(oo.dr,'restrict_columns')
-        dr.restrict_columns = oo.dr.restrict_columns;
+    if isfield(oo_.dr,'restrict_columns')
+        dr.restrict_columns = oo_.dr.restrict_columns;
     end
-    if isfield(oo.dr,'obs_var')
-        dr.obs_var = oo.dr.obs_var;
+    if isfield(oo_.dr,'obs_var')
+        dr.obs_var = oo_.dr.obs_var;
     end
 end
 
-if M.exo_nbr == 0
-    oo.exo_steady_state = [] ;
+if M_.exo_nbr == 0
+    oo_.exo_steady_state = [] ;
 end
 
-[dr.ys,M.params,info] = evaluate_steady_state(oo.steady_state,[oo.exo_steady_state; oo.exo_det_steady_state],M,options,~options.steadystate.nocheck);
+[dr.ys,M_.params,info] = evaluate_steady_state(oo_.steady_state,[oo_.exo_steady_state; oo_.exo_det_steady_state],M_,options_,~options_.steadystate.nocheck);
 
 if info(1)
-    oo.dr = dr;
+    oo_.dr = dr;
     return
 end
 
-if options.loglinear
+if options_.loglinear
     threshold = 1e-16;
     % Find variables with non positive steady state. Skip auxiliary
     % variables for lagges/leaded exogenous variables
-    idx = find(dr.ys(get_all_variables_but_lagged_leaded_exogenous(M))<threshold);
+    idx = find(dr.ys(get_all_variables_but_lagged_leaded_exogenous(M_))<threshold);
     if ~isempty(idx)
-        if options.debug
-            variables_with_non_positive_steady_state = M.endo_names{idx};
+        if options_.debug
+            variables_with_non_positive_steady_state = M_.endo_names{idx};
             skipline()
             fprintf('You are attempting to simulate/estimate a loglinear approximation of a model, but\n')
             fprintf('the steady state level of the following variables is not strictly positive:\n')
@@ -111,5 +111,5 @@ if options.loglinear
     end
 end
 
-[dr, info] = stochastic_solvers(dr, check_flag, M, options, oo);
-oo.dr = dr;
+[dr, info] = stochastic_solvers(dr, check_flag, M_, options_, oo_);
+oo_.dr = dr;
diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m
index 742b19f7ed5bdcd2108bd85049f0223e6a53fc91..864c0bd789d0e595b7cf3f3923eb263377be39e2 100644
--- a/matlab/stoch_simul.m
+++ b/matlab/stoch_simul.m
@@ -1,6 +1,20 @@
 function [info, oo_, options_, M_] = stoch_simul(M_, options_, oo_, var_list)
+%[info, oo_, options_, M_] = stoch_simul(M_, options_, oo_, var_list)
+% Computes the stochastic simulations
+%
+% INPUTS
+% - M_            [structure]       Matlab's structure describing the model
+% - options_      [structure]       Matlab's structure describing the current options
+% - oo_           [structure]       Matlab's structure containing the results
+% - var_list     [character array]  list of endogenous variables specified
+%
+% OUTPUTS
+% - info          [integer]       scalar or vector, error code.
+% - oo            [structure]     Matlab's structure containing the results (oo_).
+% - options_      [structure]       Matlab's structure describing the current options
+% - M             [structure]     Matlab's structure describing the model (M_).
 
-% Copyright © 2001-2021 Dynare Team
+% Copyright © 2001-2023 Dynare Team
 %
 % This file is part of Dynare.
 %