diff --git a/matlab/k_order_pert.m b/matlab/k_order_pert.m
index e140fd6920e35e9896a83e3bf1e49248fe274f67..02ad4bb6bde062aab7921852ba057e75ce3fc1b8 100644
--- a/matlab/k_order_pert.m
+++ b/matlab/k_order_pert.m
@@ -43,6 +43,10 @@ for i = 0:order
     gname = [ 'g_' num2str(i) ];
     dr.(gname) = dynpp_derivs.(gname);
 end
+if options.pruning
+   dr.pruning = dynpp_derivs.pruning;
+end
+
 
 % Fill equivalent Dynare matrices (up to 3rd order)
 
diff --git a/matlab/simult_.m b/matlab/simult_.m
index 37ac55d2e17bf18c31e0cf137f9a54c3d9fc9ddd..ed7f2c096fabb7dc60fa924ce64a1967c7e4f9e7 100644
--- a/matlab/simult_.m
+++ b/matlab/simult_.m
@@ -52,16 +52,15 @@ if ~options_.k_order_solver || (options_.k_order_solver && options_.pruning) %if
     end
 end
 
-if options_.k_order_solver && ~options_.pruning % Call dynare++ routines.
+if options_.k_order_solver
     if options_.order~=iorder
         error(['The k_order_solver requires the specified approximation order to be '...
                 'consistent with the one used for computing the decision rules'])
     end
-    y_start=y_(:,1); %store first period required for output
     y_ = k_order_simul(iorder,M_.nstatic,M_.npred,M_.nboth,M_.nfwrd,exo_nbr, ...
-                       y_start(dr.order_var,:),ex_',dr.ys(dr.order_var),dr);
+                       y0(dr.order_var,:),ex_',dr.ys(dr.order_var),dr, ...
+                       options_.pruning);
     y_(dr.order_var,:) = y_;
-    y_=[y_start y_];
 else
     k2 = dr.kstate(find(dr.kstate(:,2) <= M_.maximum_lag+1),[1 2]);
     k2 = k2(:,1)+(M_.maximum_lag+1-k2(:,2))*endo_nbr;
@@ -166,6 +165,6 @@ else
             yhat3 = yhat3(ipred);
         end
       otherwise
-          error(['pruning not available for order = ' int2str(iorder)])
+          error(['k_order_solver required for pruning at order = ' int2str(iorder)])
     end
 end
diff --git a/mex/build/libkordersim.am b/mex/build/libkordersim.am
index 1199cce2eb50c3ff765c88076b67e0736991dff1..f3c53ddc7ab3e918255ef5eae8460b0dd0f36a54 100644
--- a/mex/build/libkordersim.am
+++ b/mex/build/libkordersim.am
@@ -6,6 +6,7 @@ nodist_libkordersim_a_SOURCES = \
 	pascal.f08 \
 	sort.f08 \
 	partitions.f08 \
+	tensors.f08 \
 	simulation.f08 \
 	struct.f08 \
 	pthread.F08
@@ -24,7 +25,10 @@ sort.mod: sort.o
 partitions.o: sort.mod pascal.mod
 partitions.mod: partitions.o
 
-simulation.o: partitions.mod lapack.mod
+tensors.o: partitions.mod
+tensors.mod: tensors.o
+
+simulation.o: tensors.mod lapack.mod matlab_mex.mod
 simulation.mod: simulation.o
 
 %.f08: $(top_srcdir)/../../sources/libkordersim/%.f08
diff --git a/mex/sources/folded_to_unfolded_dr/mexFunction.f08 b/mex/sources/folded_to_unfolded_dr/mexFunction.f08
index 988c6c2dc63c4392a704fb1d92048cd2a2a5044a..50b5eb90458ee63caa76c9d012910d7162543dd6 100644
--- a/mex/sources/folded_to_unfolded_dr/mexFunction.f08
+++ b/mex/sources/folded_to_unfolded_dr/mexFunction.f08
@@ -27,8 +27,8 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
    type(c_ptr), dimension(*), intent(in), target :: prhs
    type(c_ptr), dimension(*), intent(out) :: plhs
    integer(c_int), intent(in), value :: nlhs, nrhs
-   type(c_ptr) :: M_mx, options_mx, dr_mx, tmp, g
-   type(pol), dimension(:), allocatable, target :: fdr
+   type(c_ptr) :: M_mx, options_mx, dr_mx, pruning_mx, tmp, g
+   type(tensor), dimension(:), allocatable, target :: fdr
    integer :: order, npred, nboth, nstatic, nfwrd, endo_nbr, exo_nbr, nys, nvar 
    type(pascal_triangle) :: p
    type(uf_matching), dimension(:), allocatable :: matching 
@@ -72,16 +72,16 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
       end if
       m = int(mxGetM(tmp))
       n = int(mxGetN(tmp))
-      allocate(fdr(d)%g(m,n))
-      fdr(d)%g = reshape(mxGetPr(tmp), [m,n])
+      allocate(fdr(d)%m(m,n))
+      fdr(d)%m = reshape(mxGetPr(tmp), [m,n])
    end do
 
    plhs(1) = mxCreateStructMatrix(1_mwSize, 1_mwSize, fieldnames)
    g = mxCreateDoubleMatrix(int(endo_nbr, mwSize), 1_mwSize, mxREAL)
-   mxGetPr(g) = reshape(fdr(0)%g, [size(fdr(0)%g)])
+   mxGetPr(g) = reshape(fdr(0)%m, [size(fdr(0)%m)])
    call mxSetField(plhs(1), 1_mwIndex, "g_0", g)
    g = mxCreateDoubleMatrix(int(endo_nbr, mwSize), int(nvar, mwSize), mxREAL)
-   mxGetPr(g) = reshape(fdr(1)%g, [size(fdr(1)%g)])
+   mxGetPr(g) = reshape(fdr(1)%m, [size(fdr(1)%m)])
    call mxSetField(plhs(1), 1_mwIndex, "g_1", g)
 
    if (order > 1) then
@@ -93,7 +93,7 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
          allocate(matching(d)%folded(nvar**d))
          call fill_folded_indices(matching(d)%folded, nvar, d, p) 
          g = mxCreateDoubleMatrix(int(endo_nbr, mwSize), int(nvar**d, mwSize), mxREAL)
-         mxGetPr(g) = reshape(fdr(d)%g(:,matching(d)%folded), [size(fdr(d)%g(:,matching(d)%folded))])
+         mxGetPr(g) = reshape(fdr(d)%m(:,matching(d)%folded), [size(fdr(d)%m(:,matching(d)%folded))])
          call mxSetField(plhs(1), 1_mwIndex, trim(fieldnames(d)), g)
       end do
    end if
diff --git a/mex/sources/k_order_mean/mexFunction.f08 b/mex/sources/k_order_mean/mexFunction.f08
index d34c7fe490d74a39ba4714ec17b6747c7a67dceb..6fd41d245a38ea1fd11dd6aaf30c664b4eeab227 100644
--- a/mex/sources/k_order_mean/mexFunction.f08
+++ b/mex/sources/k_order_mean/mexFunction.f08
@@ -46,13 +46,13 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
    integer(c_int), intent(in), value :: nlhs, nrhs
    type(c_ptr) :: order_mx, nstatic_mx, npred_mx, nboth_mx, nfwrd_mx, nexog_mx, order_moment_mx, nburn_mx, &
    yhat_start_mx, shocks_mx, ysteady_mx, dr_mx, tmp
-   type(pol), dimension(:), allocatable, target :: fdr, udr
+   type(tensor), dimension(:), allocatable, target :: fdr, udr
    integer :: order, nstatic, npred, nboth, nfwrd, exo_nbr, endo_nbr, nys, nvar, nper, nburn, order_moment
    real(real64), dimension(:,:), allocatable :: shocks, sim
    real(real64), dimension(:), allocatable :: dyu, mean
    real(real64), dimension(:), pointer, contiguous :: ysteady, yhat_start
    type(pascal_triangle) :: p
-   type(horner), dimension(:), allocatable :: h
+   type(tensor), dimension(:), allocatable :: h
    integer :: i, t, d, m, n
    character(kind=c_char, len=10) :: fieldname
 
@@ -155,12 +155,12 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
       end if
       m = int(mxGetM(tmp))
       n = int(mxGetN(tmp))
-      allocate(fdr(i)%g(m,n), udr(i)%g(endo_nbr, nvar**i), h(i)%c(endo_nbr, nvar**i))
-      fdr(i)%g(1:m,1:n) = reshape(mxGetPr(tmp), [m,n])
+      allocate(fdr(i)%m(m,n), udr(i)%m(endo_nbr, nvar**i), h(i)%m(endo_nbr, nvar**i))
+      fdr(i)%m(1:m,1:n) = reshape(mxGetPr(tmp), [m,n])
    end do
 
-   udr(0)%g = fdr(0)%g
-   udr(1)%g = fdr(1)%g
+   udr(0)%m = fdr(0)%m
+   udr(1)%m = fdr(1)%m
    if (order > 1) then
       ! Compute the useful binomial coefficients from Pascal's triangle
       p = pascal_triangle(nvar+order-1)
@@ -170,7 +170,7 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
         do d=2,order
            allocate(matching(d)%folded(nvar**d))
            call fill_folded_indices(matching(d)%folded, nvar, d, p)
-           udr(d)%g = fdr(d)%g(:,matching(d)%folded)
+           udr(d)%m = fdr(d)%m(:,matching(d)%folded)
         end do
       end block
    end if
@@ -181,15 +181,15 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
    dyu(nys+1:) = shocks(:,1) 
    ! Using the Horner algorithm to evaluate the decision rule at the chosen dyu
    call eval(h, dyu, udr, endo_nbr, nvar, order)
-   sim(:,1) = h(0)%c(:,1) + ysteady
+   sim(:,1) = h(0)%m(:,1) + ysteady
    mean = 0.
 
    ! Carrying out the simulation
    do t=2,nper
-      dyu(1:nys) = h(0)%c(nstatic+1:nstatic+nys,1) 
+      dyu(1:nys) = h(0)%m(nstatic+1:nstatic+nys,1) 
       dyu(nys+1:) = shocks(:,t)
       call eval(h, dyu, udr, endo_nbr, nvar, order)
-      sim(:,t) = h(0)%c(:,1) + ysteady
+      sim(:,t) = h(0)%m(:,1) + ysteady
       if (t > nburn) then
          mean = mean + sim(:,t)**order_moment
       end if
diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cc b/mex/sources/k_order_perturbation/k_order_perturbation.cc
index 8dd06af72e5d5b14805d69cfe94d5ef9d519096e..68f194164822b951ba71706e1c9e854aaba1523e 100644
--- a/mex/sources/k_order_perturbation/k_order_perturbation.cc
+++ b/mex/sources/k_order_perturbation/k_order_perturbation.cc
@@ -126,6 +126,11 @@ extern "C" {
       mexErrMsgTxt("options_.debug should be a logical scalar");
     bool debug = static_cast<bool>(mxGetScalar(debug_mx));
 
+    const mxArray *pruning_mx = mxGetField(options_mx, 0, "pruning");
+    if (!(pruning_mx && mxIsLogicalScalar(pruning_mx)))
+      mexErrMsgTxt("options_.pruning should be a logical scalar");
+    bool pruning = static_cast<bool>(mxGetScalar(pruning_mx));
+
     // Extract various fields from M_
     const mxArray *fname_mx = mxGetField(M_mx, 0, "fname");
     if (!(fname_mx && mxIsChar(fname_mx) && mxGetM(fname_mx) == 1))
@@ -233,7 +238,7 @@ extern "C" {
                            dr_order, llincidence);
 
         // construct main K-order approximation class
-        Approximation app(dynare, journal, nSteps, false, qz_criterium);
+        Approximation app(dynare, journal, nSteps, false, pruning, qz_criterium);
         // run stochastic steady
         app.walkStochSteady();
 
@@ -246,7 +251,18 @@ extern "C" {
         const char *g_fieldnames_c[kOrder+1];
         for (int i = 0; i <= kOrder; i++)
           g_fieldnames_c[i] = g_fieldnames[i].c_str();
-        plhs[0] = mxCreateStructMatrix(1, 1, kOrder+1, g_fieldnames_c);
+        
+        if (pruning)
+          {
+            std::vector<std::string> g_fieldnames_pruning(g_fieldnames);   
+            g_fieldnames_pruning.emplace_back("pruning");
+            const char *g_fieldnames_pruning_c[kOrder+2];
+            std::copy_n(g_fieldnames_c, kOrder+1, g_fieldnames_pruning_c);
+            g_fieldnames_pruning_c[kOrder+1] = g_fieldnames_pruning.back().c_str();
+            plhs[0] = mxCreateStructMatrix(1, 1, kOrder+2, g_fieldnames_pruning_c);
+          }
+        else
+            plhs[0] = mxCreateStructMatrix(1, 1, kOrder+1, g_fieldnames_c);
 
         // Fill that structure
         for (int i = 0; i <= kOrder; i++)
@@ -256,7 +272,27 @@ extern "C" {
             const ConstVector &vec = t.getData();
             assert(vec.skip() == 1);
             std::copy_n(vec.base(), vec.length(), mxGetPr(tmp));
-            mxSetField(plhs[0], 0, ("g_" + std::to_string(i)).c_str(), tmp);
+            mxSetField(plhs[0], 0, g_fieldnames_c[i], tmp);
+          }
+
+        // Filling the output elements for pruning 
+        if (pruning)
+          {
+            const UnfoldDecisionRule &udr_pruning = app.getUnfoldDecisionRulePruning();
+
+            mxArray *dr_pruning = mxCreateStructMatrix(1, 1, kOrder+1, g_fieldnames_c);
+            mxSetField(plhs[0], 0, "pruning", dr_pruning);
+
+            // Fill that structure
+            for (int i = 0; i <= kOrder; i++)
+              {
+                const UFSTensor &t = udr_pruning.get(Symmetry{i});
+                mxArray *tmp = mxCreateDoubleMatrix(t.nrows(), t.ncols(), mxREAL);
+                const ConstVector &vec = t.getData();
+                assert(vec.skip() == 1);
+                std::copy_n(vec.base(), vec.length(), mxGetPr(tmp));
+                mxSetField(dr_pruning, 0, g_fieldnames_c[i], tmp);
+              }
           }
 
         if (nlhs > 1)
diff --git a/mex/sources/k_order_simul/mexFunction.f08 b/mex/sources/k_order_simul/mexFunction.f08
index 30dc61d711d327b39e741c58afbb05eac95ba68c..8778fface920de394f87f4a3997f2bf979e47644 100644
--- a/mex/sources/k_order_simul/mexFunction.f08
+++ b/mex/sources/k_order_simul/mexFunction.f08
@@ -26,31 +26,26 @@
 !       shocks   matrix of shocks (nexog x number of period)
 !       ysteady  full vector of decision rule's steady
 !       dr       structure containing matrices of derivatives (g_0, g_1,…)
+!       pruning  boolean stating whether the simulation should be pruned
 !  output:
 !       res      simulated results
 
 subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
-   use iso_fortran_env
    use iso_c_binding
    use struct
-   use matlab_mex
-   use partitions
    use simulation
    implicit none (type, external)
 
    type(c_ptr), dimension(*), intent(in), target :: prhs
    type(c_ptr), dimension(*), intent(out) :: plhs
    integer(c_int), intent(in), value :: nlhs, nrhs
-   type(c_ptr) :: order_mx, nstatic_mx, npred_mx, nboth_mx, nfwrd_mx, nexog_mx, ystart_mx, shocks_mx, ysteady_mx, dr_mx, tmp
-   type(pol), dimension(:), allocatable, target :: fdr, udr
+   type(c_ptr) :: order_mx, nstatic_mx, npred_mx, nboth_mx, nfwrd_mx, &
+                 &nexog_mx, ystart_mx, shocks_mx, ysteady_mx, dr_mx, &
+                 &pruning_mx 
    integer :: order, nstatic, npred, nboth, nfwrd, exo_nbr, endo_nbr, nys, nvar, nper
-   real(real64), dimension(:,:), allocatable :: shocks, sim
-   real(real64), dimension(:), allocatable :: ysteady_pred, ystart_pred, dyu
-   real(real64), dimension(:), pointer, contiguous :: ysteady, ystart
-   type(pascal_triangle) :: p
-   type(horner), dimension(:), allocatable :: h
-   integer :: i, t, d, m, n
-   character(kind=c_char, len=10) :: fieldname
+   real(real64), dimension(:), allocatable :: dy
+   real(real64), pointer, contiguous :: ysteady(:), ystart(:), sim(:,:), shocks(:,:)
+   logical :: pruning
 
    order_mx = prhs(1)
    nstatic_mx = prhs(2)
@@ -62,10 +57,11 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
    shocks_mx = prhs(8)
    ysteady_mx = prhs(9)
    dr_mx = prhs(10)
+   pruning_mx = prhs(11)
 
    ! Checking the consistence and validity of input arguments
-   if (nrhs /= 10 .or. nlhs /= 1) then
-      call mexErrMsgTxt("Must have exactly 10 inputs and 1 output")
+   if (nrhs /= 11 .or. nlhs /= 1) then
+      call mexErrMsgTxt("Must have exactly 11 inputs and 1 output")
    end if
    if (.not. (mxIsScalar(order_mx)) .and. mxIsNumeric(order_mx)) then
       call mexErrMsgTxt("1st argument (order) should be a numeric scalar")
@@ -99,6 +95,9 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
    if (.not. mxIsStruct(dr_mx)) then
       call mexErrMsgTxt("10th argument (dr) should be a struct")
    end if
+   if (.not. (mxIsLogicalScalar(pruning_mx))) then
+      call mexErrMsgTxt("11th argument (pruning) should be a logical scalar")
+   end if
 
    ! Converting inputs in Fortran format
    order = int(mxGetScalar(order_mx))
@@ -109,6 +108,7 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
    exo_nbr = int(mxGetScalar(nexog_mx))
    endo_nbr = nstatic+npred+nboth+nfwrd
    nys = npred+nboth
+   pruning = mxGetScalar(pruning_mx) == 1._c_double
    nvar = nys+exo_nbr
 
    if (endo_nbr /= int(mxGetM(ystart_mx))) then
@@ -120,63 +120,32 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
       call mexErrMsgTxt("shocks should have nexog rows")
    end if
    nper = int(mxGetN(shocks_mx))
-   allocate(shocks(exo_nbr,nper))
-   shocks = reshape(mxGetPr(shocks_mx),[exo_nbr,nper])
+   shocks(1:exo_nbr,1:nper) => mxGetPr(shocks_mx)
 
    if (.not. (int(mxGetM(ysteady_mx)) == endo_nbr)) then
       call mexErrMsgTxt("ysteady should have nstat+npred+nboth+nforw rows")
    end if
    ysteady => mxGetPr(ysteady_mx)
+   ! Initial value for between the states' starting value and the states' 
+   ! steady-state value
+   dy = ystart(nstatic+1:nstatic+nys)-ysteady(nstatic+1:nstatic+nys)
 
-   allocate(h(0:order), fdr(0:order), udr(0:order)) 
-   do i = 0, order
-      write (fieldname, '(a2, i1)') "g_", i
-      tmp = mxGetField(dr_mx, 1_mwIndex, trim(fieldname))
-      if (.not. (c_associated(tmp) .and. mxIsDouble(tmp) .and. .not. mxIsComplex(tmp) .and. .not. mxIsSparse(tmp))) then
-         call mexErrMsgTxt(trim(fieldname)//" is not allocated in dr")
+   if (pruning) then
+      dr_mx = mxGetField(dr_mx, 1_mwIndex, "pruning")
+      if (.not. mxIsStruct(dr_mx)) then
+         call mexErrMsgTxt("dr.pruning should be a struct")
       end if
-      m = int(mxGetM(tmp))
-      n = int(mxGetN(tmp))
-      allocate(fdr(i)%g(m,n), udr(i)%g(endo_nbr, nvar**i), h(i)%c(endo_nbr, nvar**i))
-      fdr(i)%g(1:m,1:n) = reshape(mxGetPr(tmp), [m,n])
-   end do
-
-   udr(0)%g = fdr(0)%g
-   udr(1)%g = fdr(1)%g
-   if (order > 1) then
-      ! Compute the useful binomial coefficients from Pascal's triangle
-      p = pascal_triangle(nvar+order-1)
-      block
-        type(uf_matching), dimension(2:order) :: matching
-        ! Pinpointing the corresponding offsets between folded and unfolded tensors
-        do d=2,order
-           allocate(matching(d)%folded(nvar**d))
-           call fill_folded_indices(matching(d)%folded, nvar, d, p)
-           udr(d)%g = fdr(d)%g(:,matching(d)%folded)
-        end do
-      end block
    end if
 
-   allocate(dyu(nvar), ystart_pred(nys), ysteady_pred(nys), sim(endo_nbr,nper))
-   ! Getting the predetermined part of the endogenous variable vector 
-   ystart_pred = ystart(nstatic+1:nstatic+nys)
-   ysteady_pred = ysteady(nstatic+1:nstatic+nys) 
-   dyu(1:nys) = ystart_pred - ysteady_pred 
-   dyu(nys+1:) = shocks(:,1) 
-   ! Using the Horner algorithm to evaluate the decision rule at the chosen dyu
-   call eval(h, dyu, udr, endo_nbr, nvar, order)
-   sim(:,1) = h(0)%c(:,1) + ysteady
-
-   ! Carrying out the simulation
-   do t=2,nper
-      dyu(1:nys) = h(0)%c(nstatic+1:nstatic+nys,1) 
-      dyu(nys+1:) = shocks(:,t)
-      call eval(h, dyu, udr, endo_nbr, nvar, order)
-      sim(:,t) = h(0)%c(:,1) + ysteady
-   end do
-  
    ! Generating output
-   plhs(1) = mxCreateDoubleMatrix(int(endo_nbr, mwSize), int(nper, mwSize), mxREAL)
-   mxGetPr(plhs(1)) = reshape(sim, (/size(sim)/))
+   plhs(1) = mxCreateDoubleMatrix(int(endo_nbr, mwSize), int(nper+1, mwSize), mxREAL)
+   sim(1:endo_nbr,1:(nper+1)) => mxGetPr(plhs(1))
+   sim(:,1) = ystart
+
+   if (pruning) then
+      call simulate_pruning(sim, dr_mx, ysteady, dy, shocks, order, nstatic, nvar)
+   else
+      call simulate(sim, dr_mx, ysteady, dy, shocks, order, nstatic, nvar)
+   end if
 
 end subroutine mexFunction
diff --git a/mex/sources/k_order_welfare/k_order_welfare.cc b/mex/sources/k_order_welfare/k_order_welfare.cc
index 711e4d90064a2daf34128c00fc1f3f9a3cf8b16e..259f21705438d98ec3ff2fbd8b89b58de77cb42a 100644
--- a/mex/sources/k_order_welfare/k_order_welfare.cc
+++ b/mex/sources/k_order_welfare/k_order_welfare.cc
@@ -139,6 +139,11 @@ extern "C" {
       mexErrMsgTxt("options_.debug should be a logical scalar");
     bool debug = static_cast<bool>(mxGetScalar(debug_mx));
 
+    const mxArray *pruning_mx = mxGetField(options_mx, 0, "pruning");
+    if (!(pruning_mx && mxIsLogicalScalar(pruning_mx)))
+      mexErrMsgTxt("options_.pruning should be a logical scalar");
+    bool pruning = static_cast<bool>(mxGetScalar(pruning_mx));
+
     // Extract various fields from M_
     const mxArray *fname_mx = mxGetField(M_mx, 0, "fname");
     if (!(fname_mx && mxIsChar(fname_mx) && mxGetM(fname_mx) == 1))
@@ -262,7 +267,7 @@ extern "C" {
 
 
     // construct main K-order approximation class
-    Approximation app(dynare, journal, nSteps, false, qz_criterium);
+    Approximation app(dynare, journal, nSteps, false, pruning, qz_criterium);
     // run stochastic steady
     app.walkStochSteady();
 
diff --git a/mex/sources/libkorder/kord/approximation.cc b/mex/sources/libkorder/kord/approximation.cc
index 623d3f3a8e9fcef3ba53c213e25f8665f85f686e..d07a333e5e5cf66c0debffabf8260cc7fe8a7c2c 100644
--- a/mex/sources/libkorder/kord/approximation.cc
+++ b/mex/sources/libkorder/kord/approximation.cc
@@ -49,13 +49,14 @@ ZAuxContainer::getType(int i, const Symmetry &s) const
   return itype::zero;
 }
 
-Approximation::Approximation(DynamicModel &m, Journal &j, int ns, bool dr_centr, double qz_crit)
+Approximation::Approximation(DynamicModel &m, Journal &j, int ns, bool dr_centr, bool pruned_dr, double qz_crit)
   : model(m), journal(j),
     ypart(model.nstat(), model.npred(), model.nboth(), model.nforw()),
     mom(UNormalMoments(model.order(), model.getVcov())),
     nvs{ypart.nys(), model.nexog(), model.nexog(), 1},
     steps(ns),
-    dr_centralize(dr_centr), qz_criterium(qz_crit), ss(ypart.ny(), steps+1)
+    dr_centralize(dr_centr), pruning(pruned_dr),
+    qz_criterium(qz_crit), ss(ypart.ny(), steps+1)
 {
   ss.nans();
 }
@@ -69,6 +70,16 @@ Approximation::getFoldDecisionRule() const
   return *fdr;
 }
 
+/* This just returns ‘fdr_pruning’ with a check that it is created. */
+const UnfoldDecisionRule &
+Approximation::getUnfoldDecisionRulePruning() const
+{
+  KORD_RAISE_IF(!udr_pruning,
+                "Unfolded decision rule has not been created in Approximation::getUnfoldDecisionRule");
+  return *udr_pruning;
+}
+
+
 /* This just returns ‘udr’ with a check that it is created. */
 const UnfoldDecisionRule &
 Approximation::getUnfoldDecisionRule() const
@@ -209,6 +220,15 @@ Approximation::walkStochSteady()
   udr.reset();
   fdr = std::make_unique<FoldDecisionRule>(*rule_ders, ypart, model.nexog(),
                                            model.getSteady(), 1.0-sigma_so_far);
+  if (pruning)
+    {
+      fdr_pruning = std::make_unique<FoldDecisionRule>(*rule_ders, ypart, 
+                                                     model.nexog(), 
+                                                     model.getSteady(), 
+                                                     1.0-sigma_so_far,
+                                                     pruning);
+      udr_pruning = std::make_unique<UnfoldDecisionRule>(*fdr_pruning);
+    } 
   if (steps == 0 && dr_centralize)
     {
       // centralize decision rule for zero steps
diff --git a/mex/sources/libkorder/kord/approximation.hh b/mex/sources/libkorder/kord/approximation.hh
index 42b31290193e0928a43cba6ed478ddc1423b1113..e2b2ae505bf73bb36f6f8d6ee99f69fe27b5ebab 100644
--- a/mex/sources/libkorder/kord/approximation.hh
+++ b/mex/sources/libkorder/kord/approximation.hh
@@ -125,18 +125,22 @@ class Approximation
   std::unique_ptr<FGSContainer> rule_ders_s;
   std::unique_ptr<FGSContainer> rule_ders_ss;
   std::unique_ptr<FoldDecisionRule> fdr;
+  std::unique_ptr<FoldDecisionRule> fdr_pruning;
   std::unique_ptr<UnfoldDecisionRule> udr;
+  std::unique_ptr<UnfoldDecisionRule> udr_pruning;
   const PartitionY ypart;
   const FNormalMoments mom;
   IntSequence nvs;
   int steps;
   bool dr_centralize;
+  bool pruning;
   double qz_criterium;
   TwoDMatrix ss;
 public:
-  Approximation(DynamicModel &m, Journal &j, int ns, bool dr_centr, double qz_crit);
+  Approximation(DynamicModel &m, Journal &j, int ns, bool dr_centr, bool pruning, double qz_crit);
 
   const FoldDecisionRule &getFoldDecisionRule() const;
+  const UnfoldDecisionRule &getUnfoldDecisionRulePruning() const;
   const UnfoldDecisionRule &getUnfoldDecisionRule() const;
   const TwoDMatrix &
   getSS() const
diff --git a/mex/sources/libkorder/kord/decision_rule.hh b/mex/sources/libkorder/kord/decision_rule.hh
index 9d7fb732a929e928363537aa2bbb309f5b6656d2..c2ec0d3cdb1918ffdd2335bb98325f14751d988c 100644
--- a/mex/sources/libkorder/kord/decision_rule.hh
+++ b/mex/sources/libkorder/kord/decision_rule.hh
@@ -129,6 +129,16 @@ public:
   {
     fillTensors(g, sigma);
   }
+  DecisionRuleImpl(const _Tg &g, const PartitionY &yp, int nuu,
+                   const ConstVector &ys, double sigma, bool pruning)
+    : ctraits<t>::Tpol(yp.ny(), yp.nys()+nuu), ysteady(ys), ypart(yp), nu(nuu)
+  {
+    if (pruning)
+      fillTensorsPruning(g);
+    else
+      fillTensors(g, sigma);
+  }
+
   DecisionRuleImpl(const _TW &W, int nys, int nuu,
                    const ConstVector &ys)
     : ctraits<t>::Tpol(1, nys+nuu), ysteady(ys), nu(nuu)
@@ -162,6 +172,7 @@ public:
   }
 protected:
   void fillTensors(const _Tg &g, double sigma);
+  void fillTensorsPruning(const _Tg &g);
   void fillTensors(const _TW &W, int nys);
   void centralize(const DecisionRuleImpl &dr);
 public:
@@ -247,6 +258,43 @@ DecisionRuleImpl<t>::fillTensors(const _Tg &g, double sigma)
     }
 }
 
+template<Storage t>
+void
+DecisionRuleImpl<t>::fillTensorsPruning(const _Tg &g)
+{
+  IntSequence tns{ypart.nys(), nu, 1};
+  int dfact = 1;
+  for (int d = 0; d <= g.getMaxDim(); d++, dfact *= d)
+    {
+      auto g_yusd = std::make_unique<_Ttensym>(ypart.ny(), ypart.nys()+nu+1, d);
+      g_yusd->zeros();
+      // fill tensor of ‘g_yusd’ of dimension ‘d’
+      /* 
+        Here we have to fill the tensor [g_(yuσ)ᵈ]. So we go through all pairs
+        (i,j,k) such that i+j+k=d. We weight it with 1/(i+j+k)! The factorial
+        is denoted dfact.
+      */
+      for (int i = 0; i <= d; i++)
+        {
+          for (int j = 0; j <= d-i; j++)
+            {
+              int k = d-i-j;
+              _Ttensor tmp(ypart.ny(),
+                           TensorDimens(Symmetry{i, j, k}, tns));
+              tmp.zeros();
+              if (Symmetry sym{i, j, 0, k}; g.check(sym))
+                {
+                  double mult = 1.0/dfact;
+                  // mexPrintf("Symmetry found: %d %d %d %.2f %d\n", i, j, k, mult, kfact);
+                  tmp.add(mult, g.get(sym));
+                }
+              g_yusd->addSubTensor(tmp);
+            }
+        }
+      this->insert(std::move(g_yusd));
+    }
+}
+
 template<Storage t>
 void
 DecisionRuleImpl<t>::fillTensors(const _TW &W, int nys)
@@ -396,6 +444,11 @@ public:
     : DecisionRuleImpl<Storage::fold>(g, yp, nuu, ys, sigma)
   {
   }
+  FoldDecisionRule(const ctraits<Storage::fold>::Tg &g, const PartitionY &yp, int nuu,
+                   const ConstVector &ys, double sigma, bool pruning)
+    : DecisionRuleImpl<Storage::fold>(g, yp, nuu, ys, sigma, pruning)
+  {
+  }
   FoldDecisionRule(const ctraits<Storage::fold>::TW &W, int nys, int nuu,
                    const ConstVector &ys)
     : DecisionRuleImpl<Storage::fold>(W, nys, nuu, ys)
diff --git a/mex/sources/libkordersim/partitions.f08 b/mex/sources/libkordersim/partitions.f08
index 878cea43c3382cc1db716dc8eb22e134618a8ee3..d3fc42b79de4033f233ab793f22b5e64a177a086 100644
--- a/mex/sources/libkordersim/partitions.f08
+++ b/mex/sources/libkordersim/partitions.f08
@@ -1,7 +1,3 @@
-! Provides subroutines to manipulate indexes representing elements of
-! a partition for a given integer
-! i.e. elements p = (α₁,…,αₘ) where each αᵢ ∈ { 0, ..., n-1 }
-
 ! Copyright © 2021-2023 Dynare Team
 !
 ! This file is part of Dynare.
@@ -54,6 +50,14 @@ module partitions
    type uf_matching
       type(integer), dimension(:), allocatable :: folded
    end type
+   
+
+   ! A type to contain the integer partitions up to integer n
+   ! with up to n parts
+   type partition_triangle
+      type(index), dimension(:), allocatable :: partition ! integer partitions
+      integer, dimension(:), allocatable :: count ! numbers of equivalent permuted partitions
+   end type partition_triangle
 
 contains
 
@@ -184,8 +188,12 @@ contains
          find = 0
       else
          i = 1
-         do while (i <= l .and. a(i) /= v)
-            i = i+1
+         do while (i <= l)
+            if (a(i) /= v) then
+               i = i+1
+            else
+               exit
+            end if
          end do
          if (i == l+1) then
             find = 0
@@ -295,7 +303,7 @@ contains
    ! The algorithm to count the number of equivalent unfolded indices
    ! works as follows. A folded index can be written as α = (x₁, ..., x₁, ..., xₚ, ..., xₚ)
    ! such that x₁ < x₂ < ... < xₚ. Denote kᵢ the number of coordinates equal to xᵢ.
-   ! The number of unfolded indices equivalent to α is c(α) = ⎛         d       ⎞
+   ! The number of unfolded indices equivalent to α is c(α) = ⎛         m       ⎞
    !                                                          ⎝ k₁, k₂, ..., kₚ ⎠
    ! Suppose j is the latest incremented coordinate.
    ! If αⱼ < n, then αⱼ' = αⱼ + 1, k(αⱼ) := k(αⱼ)-1, k(αⱼ') := k(αⱼ')+1.
@@ -347,6 +355,84 @@ contains
       end do
    end subroutine folded_offset_loop 
 
+   ! The following routine computes the partitions of all integers up to integer
+   ! n with number of parts up to n, where n is the size of the input
+   ! partition_triangle leading partitions. The partitions are stored in the
+   ! lexicographic order. Suppose we want to partition the
+   ! integer n using k parts, i.e. solve the problem 𝓟(k,n). Two cases arise:
+   ! (i) the partition begins with a one, i.e. is of the form (1, α₁, ...,
+   ! αₖ₋₁). In this case, (α₁, ..., αₖ₋₁) is a partition of n-1 with k-1 parts,
+   ! i.e. solves the problem 𝓟(k-1,n-1). (ii) Otherwise, if (α₁, ..., αₖ) is a
+   ! partition of integer n using k parts, then (α₁-1, ..., αₖ-1) is a partition
+   ! of integer n-k using k parts, i.e. solves the problem 𝓟(k,n-k). In other
+   ! words, solutions to the problem 𝓟(k,n) are (a) solutions to the problem
+   ! 𝓟(k-1,n-1) with a 1 appended up front, and (b) solutions to the problem
+   ! 𝓟(k,n-k) with a 1 added to all its indices. Denoting p(k,n) the cardinal
+   ! of 𝓟(k,n), it thus verifies p(k,n)=p(k-1,n-1)+p(k,n-k).
+   ! A partition of n with k parts can be written as 
+   ! α = (α₁, ..., α₁, ..., αₚ, ..., αₚ) such that α₁ < α₂ < ... < αₚ. 
+   ! Denote kᵢ the number of coordinates equal to αᵢ. The
+   ! number of partitions that are permuted version of α 
+   ! is c(α;k) = ⎛       k         ⎞.
+   !             ⎝ k₁, k₂, ..., kₚ ⎠
+   ! The partitions generated through (b) represent the same number of permuted
+   ! partitions. If α solves 𝓟(k,n-k), α'= α.+1 is such that c(α';k)=c(α;k).
+   ! As for (a), two cases arise. If α that solves 𝓟(k-1,n-1) is such that
+   ! α₁=1, then α'=[1 α] that solves 𝓟(k,n), then
+   ! c(α';k) = ⎛       k       ⎞ = ⎛       k-1+1       ⎞ = k*c(α;k-1)/(k₁+1)
+   !           ⎝ k₁', ..., kₚ' ⎠   ⎝ k₁+1, k₂, ..., kₚ ⎠
+   ! Otherwise, c(α';k)=k*c(α;k-1)
+   subroutine fill_partition_triangle(parts)
+      type(partition_triangle), dimension(:,:), intent(inout) :: parts 
+      integer :: n, k, p_km1_nm1, p_k_nmk, p_k_n, l
+      ! Initialization with n=1 
+      parts(1,1)%partition = [index(1,1)]
+      parts(1,1)%count = [1]
+      do n=2,size(parts,1)
+         ! print *, 'n:', n
+         ! 𝓟(n,n) unique solution is (1,...,1)
+         !                            n times
+         parts(n,n)%partition = [index(n,1)] 
+         parts(n,n)%count = [1]
+         ! 𝓟(1,n) unique solution is (n)
+         parts(1,n)%partition = [index(1,n)]
+         parts(1,n)%count = [1]
+         do k=2,n-1
+            ! print *, 'k:', k
+            p_km1_nm1 = size(parts(k-1,n-1)%partition)
+            if (2*k>n) then
+               p_k_nmk = 0
+            else   
+               p_k_nmk = size(parts(k,n-k)%partition)
+            end if
+            p_k_n = p_km1_nm1+p_k_nmk
+            ! print *, 'p_km1_nm1:', p_km1_nm1
+            ! print *, 'p_k_nmk:', p_k_nmk
+            ! print *, 'p_k_n:', p_k_n
+            allocate(parts(k,n)%partition(p_k_n))
+            allocate(parts(k,n)%count(p_k_n))
+            ! 1 appended up front to 𝓟(k-1,n-1) solutions 
+            do l=1,p_km1_nm1
+               ! print *, 'l:', l
+               allocate(parts(k,n)%partition(l)%coor(k))
+               parts(k,n)%partition(l)%coor(1) = 1
+               parts(k,n)%partition(l)%coor(2:k) = parts(k-1,n-1)%partition(l)%coor
+               if (parts(k-1,n-1)%partition(l)%coor(1) == 1) then
+                  parts(k,n)%count(l) = parts(k-1,n-1)%count(l)*k/(get_prefix_length(parts(k-1,n-1)%partition(l), k-1)+1)
+               else
+                  parts(k,n)%count(l) = parts(k-1,n-1)%count(l)*k
+               end if
+            end do
+            ! 1 added to all components of 𝓟(k,n-k) solutions
+            do l=1,p_k_nmk
+               ! print *, 'l:', l
+               parts(k,n)%partition(l+p_km1_nm1)%coor = parts(k,n-k)%partition(l)%coor+1
+               parts(k,n)%count(l+p_km1_nm1) = parts(k,n-k)%count(l)
+            end do
+         end do
+      end do
+   end subroutine fill_partition_triangle
+
 end module partitions
 
 ! gfortran -o partitions partitions.f08 pascal.f08 sort.f08
@@ -357,10 +443,11 @@ end module partitions
 !    implicit none (type, external)
 !    type(index) :: uidx, fidx, i1, i2
 !    integer, dimension(:), allocatable :: folded
-!    integer :: i, uj, n, d, j, nb_folded_idcs
+!    integer :: i, uj, n, d, j, nb_folded_idcs, k, l, m
 !    type(pascal_triangle) :: p
 !    type(index), dimension(:), allocatable :: list_folded_idcs
 !    integer, dimension(:), allocatable :: nbeq, off
+!    type(partition_triangle), allocatable, dimension(:,:) :: t
 
 !    ! Unfolded indices and offsets
 !    ! 0,0,0  1    1,0,0  10   2,0,0  19
@@ -425,4 +512,18 @@ end module partitions
 !    print '(i3)', (nbeq(i), i=1,nb_folded_idcs)
 !    print '(i4)', (off(i), i=1,nb_folded_idcs)
 
-! end program test
+!    ! Triangle of integer partitions
+!    d = 7
+!    allocate(t(d,d))
+!    call fill_partition_triangle(t)
+!    do n=1,d
+!       do k=1,n
+!          print *, '(k,n):', k, n
+!          do l=1, size(t(k,n)%partition)
+!             print '(10i2)', (t(k,n)%partition(l)%coor(m), m=1,k)
+!             print '(i2)', t(k,n)%count(l)
+!          end do
+!       end do
+!    end do
+
+! end program test
\ No newline at end of file
diff --git a/mex/sources/libkordersim/simulation.f08 b/mex/sources/libkordersim/simulation.f08
index 4a81f51a9031406f6074572814f57d20e382038a..ed50b644de81b60912dda8292601b1e61774037f 100644
--- a/mex/sources/libkordersim/simulation.f08
+++ b/mex/sources/libkordersim/simulation.f08
@@ -21,48 +21,40 @@
 
 module simulation
    use iso_fortran_env
-   use partitions
+   use tensors
    use lapack
+   use matlab_mex
    implicit none (type, external)
 
-   ! Used to store the folded decision rule tensors
-   type :: pol
-      real(real64), allocatable :: g(:,:)
-   end type pol
-
-   ! Used to store the different contracted tensors used in the Horner algorithm
-   ! type :: horner
-   !    real(real64), pointer, contiguous :: c(:,:), c_flat(:)
-   ! end type horner
-   type :: horner
-      real(real64), allocatable :: c(:,:)
-   end type horner
-
+   ! Used to store the simulated pruned state space
+   type pruned
+      real(real64), allocatable :: y(:,:)
+   end type pruned
 
 contains
 
    ! ! With MATMUL
    ! ! Horner evaluation of the polynomial with coefficients stored in udr at the point dyu
    ! subroutine eval(h, dyu, udr, ny, nvar, order)
-   !    type(horner), dimension(0:order), intent(inout) :: h
+   !    type(tensor), dimension(0:order), intent(inout) :: h
    !    real(real64), dimension(nvar), intent(in) :: dyu
-   !    type(pol), intent(in) :: udr(0:order)
+   !    type(tensor), intent(in) :: udr(0:order)
    !    integer, intent(in) :: ny, nvar, order
    !    integer :: d
    !    if (order == 1) then
-   !       h(1)%c = udr(1)%g
+   !       h(1)%c = udr(1)%m
    !    else
    !       call contract(h(order-1)%c, udr(order)%g, dyu, ny, nvar, order)
    !    end if
    !    do d=order-1,1,-1
    !       if (d > 1) then
-   !          h(d)%c = h(d)%c + udr(d)%g
+   !          h(d)%c = h(d)%c + udr(d)%m
    !          call contract(h(d-1)%c, h(d)%c, dyu, ny, nvar, d)
    !       else
-   !          h(d)%c = h(d)%c + udr(1)%g
+   !          h(d)%c = h(d)%c + udr(1)%m
    !       end if
    !    end do
-   !    h(0)%c(:,1) = matmul(h(1)%c, dyu) + udr(0)%g(:,1)
+   !    h(0)%c(:,1) = matmul(h(1)%c, dyu) + udr(0)%m(:,1)
    ! end subroutine eval
 
    ! ! Contracts the unfolded tensor t with respect to the vector x and stores the result in c
@@ -90,29 +82,30 @@ contains
 
    ! Horner evaluation of the polynomial with coefficients stored in udr at the point dyu
    subroutine eval(h, dyu, udr, ny, nvar, order)
-      type(horner), dimension(0:order), intent(inout) :: h
+      type(tensor), dimension(0:order), intent(inout) :: h
       real(real64), dimension(nvar), intent(in) :: dyu
-      type(pol), intent(in) :: udr(0:order)
+      type(tensor), intent(in) :: udr(0:order)
       integer, intent(in) :: ny, nvar, order
       integer :: d
       if (order == 1) then
-         h(1)%c = udr(1)%g
+         h(1)%m = udr(1)%m
       else
-         call contract(h(order-1)%c, udr(order)%g, dyu, ny, nvar, order)
+         call contract(h(order-1)%m, udr(order)%m, dyu, ny, nvar, order)
       end if
       do d=order-1,1,-1
          if (d > 1) then
-            h(d)%c = h(d)%c + udr(d)%g
-            call contract(h(d-1)%c, h(d)%c, dyu, ny, nvar, d)
+            h(d)%m = h(d)%m + udr(d)%m
+            call contract(h(d-1)%m, h(d)%m, dyu, ny, nvar, d)
          else
-            h(d)%c = h(d)%c + udr(1)%g
+            h(d)%m = h(d)%m + udr(1)%m
          end if
       end do
-      h(0)%c = udr(0)%g
-      call mult_vec(1.0_real64, h(1)%c, dyu, 1.0_real64, h(0)%c(:,1))
+      h(0)%m = udr(0)%m
+      call mult_vec(1.0_real64, h(1)%m, dyu, 1.0_real64, h(0)%m(:,1))
    end subroutine eval
 
-   ! Contracts the unfolded tensor t with respect to the vector x and stores the result in c
+   ! Contracts the unfolded tensor t with respect to the vector x and stores the
+   ! result in c.
    subroutine contract(c, t, x, nrows, nvar, d)
       integer, intent(in) :: nrows, nvar, d
       real(real64), dimension(nrows, nvar**d), intent(in) :: t
@@ -124,5 +117,171 @@ contains
       end do
    end subroutine contract
 
+   ! Simulates the model around the steady state using the pruned state space
+   subroutine simulate_pruning(sim, dr_mx, ysteady, dy, shocks, order, nstatic, nvar)
+      real(real64), dimension(:,:), contiguous, intent(inout) :: sim
+      type(c_ptr), intent(in) :: dr_mx
+      real(real64), contiguous, intent(in) :: ysteady(:), dy(:) 
+      real(real64), contiguous, target, intent(in) :: shocks(:,:)
+      integer, intent(in) :: order, nstatic, nvar
+
+      real(real64), dimension(:,:), allocatable :: phat
+      real(real64), dimension(:), contiguous, pointer :: u
+      type(tensor), dimension(:), allocatable :: psim
+      type(partition_triangle), dimension(:,:), allocatable, target :: part
+      real(real64) :: prod, sum_part
+      type(partition_triangle), pointer :: part_set
+      character(kind=c_char, len=10) :: fieldname
+      type(c_ptr) :: g_q
+      type(unfolded_tensor), dimension(:), allocatable :: g
+      integer :: q, l, j, m, r, t, c, endo_nbr, nys, nper
+      ! Import the MATLAB decision rule matrices
+      q = 1
+      allocate(g(1:order))
+      do while (q <= order)
+         write (fieldname, '(a2, i1)') "g_", q
+         g_q = mxGetField(dr_mx, 1_mwIndex, trim(fieldname))
+         if (.not. (c_associated(g_q) .and. mxIsDouble(g_q) .and. .not. mxIsComplex(g_q) .and. .not. mxIsSparse(g_q))) then
+            call mexErrMsgTxt(trim(fieldname)//" is not allocated in dr.pruning")
+         end if
+         g(q)%m(1:mxGetM(g_q),1:mxGetN(g_q)) => mxGetPr(g_q)
+         q = q+1
+      end do
+      ! Initialize useful variables
+      nys = size(dy)
+      endo_nbr = size(ysteady)
+      nper = size(sim,2)
+      allocate(psim(order), phat(nys,order), part(order,order))
+      allocate(psim(1)%m(endo_nbr, 1))
+      phat(:,1) = dy
+      do l=2,order
+         allocate(psim(l)%m(endo_nbr, 1))
+         psim(l)%m(:,1) = 0._real64
+         phat(:,l) = 0._real64
+      end do
+      call fill_list_unfolded_tensor(g, nys, nvar+1)
+      call fill_partition_triangle(part)
+      ! Loop over periods
+      do t=2,nper
+         u => shocks(:,t-1)
+         do l=1,order
+            psim(l)%m(:,1) = 0._real64
+         end do
+         ! Go over the [gᵥˡ] folded tensors
+         do l=1,order
+            ! Go over the offsets of the matrix spanning [gᵥˡ]
+            do j=1,size(g(l)%ind)
+               if (g(l)%s(j) > 0) then
+                  ! Compute the contribution of g_l-column-multiplied
+                  ! terms to x_q when the index of the g_l columns
+                  ! involves at least one state variable
+                  ! Contributions are non-zero if q ≥ l
+                  do q=l,order
+                     ! 1. Get the sum-over-integer-partition terms
+                     sum_part = 0._real64
+                     part_set => part(g(l)%s(j),q-(l-g(l)%s(j)))
+                     do r=1,size(part_set%partition)
+                        prod = real(part_set%count(r),real64)
+                        c = 1
+                        do m=1,l
+                           if (g(l)%ind(j)%coor(m) <= nys) then
+                              prod = prod*phat(g(l)%ind(j)%coor(m),&
+                              &part_set%partition(r)%coor(c))
+                              c = c+1
+                           end if
+                           if ((g(l)%ind(j)%coor(m) > nys) .and. (g(l)%ind(j)%coor(m) <= nvar)) then
+                              prod = prod*u(g(l)%ind(j)%coor(m)-nys)
+                           end if
+                        end do
+                        sum_part = sum_part+prod
+                     end do
+                     ! 2. Multiply the sum-over-integer-partition terms
+                     ! by the considered column [gᵥˡ] 
+                     psim(q)%m(:,1) = psim(q)%m(:,1)+g(l)%m(:,j)*sum_part
+                  end do
+               else
+                  ! Compute the contribution of g_l-column-multiplied
+                  ! terms to x_q when the index of the g_l columns
+                  ! involves no state variables. Such a contribution only 
+                  ! exists when q=l.
+                  prod = 1._real64
+                  do m=1,l
+                     if ((g(l)%ind(j)%coor(m) > nys) .and. (g(l)%ind(j)%coor(m) <= nvar)) then
+                        prod = prod*u(g(l)%ind(j)%coor(m)-nys)
+                     end if
+                  end do
+                  psim(l)%m(:,1) = psim(l)%m(:,1)+g(l)%m(:,j)*prod
+               end if
+            end do
+         end do
+         ! Add l-order terms to the output simulation for period t
+         ! for 1 ≤ l ≤ order
+         sim(:,t) = ysteady
+         ! Prepare next iteration
+         do l=1,order
+            phat(:,l) = psim(l)%m(nstatic+1:nstatic+nys,1)
+            sim(:,t) = sim(:,t)+psim(l)%m(:,1)
+         end do
+      end do
+   end subroutine simulate_pruning
+
+   ! Simulates the a model around the steady state
+   subroutine simulate(sim, dr_mx, ysteady, dy, shocks, order, nstatic, nvar)
+      real(real64), dimension(:,:), contiguous, intent(inout) :: sim
+      type(c_ptr), intent(in) :: dr_mx
+      real(real64), contiguous, intent(in) :: ysteady(:), dy(:), shocks(:,:)
+      integer, intent(in) :: order, nstatic, nvar
+      character(kind=c_char, len=10) :: fieldname
+      type(c_ptr) :: g_d
+      type(tensor), dimension(:), allocatable :: fg, ug, h
+      integer :: d, endo_nbr, nys, t
+      type(uf_matching), dimension(:), allocatable :: matching
+      real(real64), dimension(:), allocatable :: dyu
+      type(pascal_triangle) :: p
+      ! Import the MATLAB decision rule matrices
+      d = 0
+      allocate(fg(0:order))
+      do while (d <= order)
+         write (fieldname, '(a2, i1)') "g_", d
+         g_d = mxGetField(dr_mx, 1_mwIndex, trim(fieldname))
+         if (.not. (c_associated(g_d) .and. mxIsDouble(g_d) .and. .not. mxIsComplex(g_d) .and. .not. mxIsSparse(g_d))) then
+            call mexErrMsgTxt(trim(fieldname)//" is not allocated in dr")
+         end if
+         fg(d)%m(1:mxGetM(g_d),1:mxGetN(g_d)) => mxGetPr(g_d)
+         d = d+1
+      end do
+      ! Put decision rules matrices in the unfolded form
+      ! Pinpointing the corresponding offsets between folded and unfolded
+      ! tensors
+      allocate(h(0:order), ug(0:order), matching(2:order))
+      endo_nbr = size(ysteady)
+      do d=0,1
+         allocate(h(d)%m(endo_nbr,nvar**d))
+         ug(d)%m => fg(d)%m
+      end do
+      ! Compute the useful binomial coefficients from Pascal's triangle
+      p = pascal_triangle(nvar+order-1)
+      do d=2,order
+         allocate(h(d)%m(endo_nbr,nvar**d), ug(d)%m(endo_nbr, nvar**d), matching(d)%folded(nvar**d))
+         call fill_folded_indices(matching(d)%folded, nvar, d, p)
+         ug(d)%m = fg(d)%m(:,matching(d)%folded)
+      end do
+      allocate(dyu(nvar))
+      ! Getting the predetermined part of the endogenous variable vector 
+      nys = size(dy)
+      dyu(1:nys) = dy 
+      ! Carrying out the simulation
+      do t=2,size(sim,2)
+         dyu(nys+1:) = shocks(:,t-1)
+         ! Using the Horner algorithm to evaluate the decision rule at the 
+         ! chosen dyu
+         call eval(h, dyu, ug, endo_nbr, nvar, order)
+         sim(:,t) = h(0)%m(:,1) + ysteady
+         dyu(1:nys) = h(0)%m(nstatic+1:nstatic+nys,1)
+      end do
+   end subroutine simulate
+
+
+
 end module simulation
 
diff --git a/mex/sources/libkordersim/tensors.f08 b/mex/sources/libkordersim/tensors.f08
new file mode 100644
index 0000000000000000000000000000000000000000..728c624fb95dc227fec64d25e9b27c5b1c337936
--- /dev/null
+++ b/mex/sources/libkordersim/tensors.f08
@@ -0,0 +1,304 @@
+! Copyright © 2021-2023 Dynare Team
+!
+! This file is part of Dynare.
+!
+! Dynare is free software: you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! Dynare is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License
+! along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+
+module tensors
+   use iso_c_binding
+   use partitions
+   ! use matlab_mat
+   implicit none (type, external)
+
+   ! A type to contain a folded or unfolded tensor [gᵥᵐ]
+   type tensor
+      real(real64), pointer, contiguous, dimension(:,:) :: m
+      ! real(real64), dimension(:,:), allocatable :: m
+   end type tensor 
+
+   ! A type to contain the unfolded tensor [gᵥᵐ]. For each offset j, ind(j)
+   ! contains the associated unfolded index and s(j) stores the number of state 
+   ! variables in the index, i.e. if v = (x,u[,σ]) where x is the s-sized state 
+   ! variables vector and ind(j) = (α₁,…,αₘ), s(j) stores the number of 
+   ! coordinates αᵢ such that  αᵢ ≤ s.
+   type unfolded_tensor
+      type(index), dimension(:), allocatable :: ind
+      integer, dimension(:), allocatable :: s
+      real(real64), dimension(:,:), pointer, contiguous :: m
+   end type unfolded_tensor
+
+   ! A type to contain the folded tensor [gᵥᵐ]. For each offset j, ind(j)
+   ! contains the associated folded index, count(j) contains the number of
+   ! equivalent unfolded indices and s(j) stores the number of state variables
+   ! in the index, i.e. if v = (x,u[,σ]) where x is the s-sized state variables
+   ! vector and ind(j) = (α₁,…,αₘ), s(j) stores the number of coordinates αᵢ
+   ! such that  αᵢ ≤ s.
+   type folded_tensor
+      type(index), dimension(:), allocatable :: ind
+      integer, dimension(:), allocatable :: count
+      integer, dimension(:), allocatable :: s
+      real(real64), dimension(:,:), pointer, contiguous :: m
+   end type folded_tensor
+
+contains
+
+   ! Fills the ind and s members of a list of unfolded tensors g
+   ! with a number of state variables ns. More precisely, g contains the
+   ! information associated with the unfolded tensors [gᵥᵐ] with m=1,...,q
+   ! 1. Filling g%ind
+   ! The unfolded indices of [gᵥᵐ⁺¹] are the unfolded indices of [gᵥᵐ]
+   ! with all numbers 1,...,n appended upfront. In other words, any
+   ! index of [gᵥᵐ⁺¹] writes (α₁,α) with α₁∈ { 1, ..., n } and α 
+   ! an index of [gᵥᵐ]
+   ! 2. Filling g%s
+   ! s(α₁,α) is s(α)+1 if α₁ ≤ ns and s(α) otherwise 
+   subroutine fill_list_unfolded_tensor(g, s, n)
+      type(unfolded_tensor), dimension(:), intent(inout) :: g
+      integer, intent(in) :: s, n ! number of state variables and size of v
+      integer :: q, m, j, k, l, length
+      q = size(g)
+      ! Initialization
+      m = 1
+      length = n
+      allocate(g(1)%ind(length), g(1)%s(length))
+      do j=1,n
+         g(1)%ind(j) = index(1,j)
+         if (j<=s) then
+            g(1)%s(j) = 1
+         else
+            g(1)%s(j) = 0
+         end if
+      end do
+      ! Transmission
+      do m=2,q
+         length = n*length
+         allocate(g(m)%ind(length), g(m)%s(length))
+         l = 1
+         do j=1,n
+            do k=1,size(g(m-1)%ind)
+               g(m)%ind(l) = index(m)
+               g(m)%ind(l)%coor(1) = j
+               g(m)%ind(l)%coor(2:m) = g(m-1)%ind(k)%coor
+               if (j<=s) then
+                  g(m)%s(l) = g(m-1)%s(k)+1
+               else
+                  g(m)%s(l) = g(m-1)%s(k)
+               end if
+               l = l+1
+            end do
+         end do
+      end do
+   end subroutine fill_list_unfolded_tensor
+
+   ! Allocates the members ind, count and s of a tensor [gᵥᵐ] denoted g 
+   ! using the Pascal triangle p. n is the size of v.
+   subroutine allocate_folded_tensor(g, m, n, p)
+      type(folded_tensor), intent(inout) :: g
+      integer, intent(in) :: m, n ! 
+      type(pascal_triangle),intent(in) :: p
+      integer :: d
+      d = get(m,n+m-1,p)
+      allocate(g%ind(d), g%count(d), g%s(d))
+   end subroutine allocate_folded_tensor
+   
+   ! Fills the already allocated members ind, count and s of a tensor [gᵥᵐ] 
+   ! denoted g using the Pascal triangle p. n is the size of v.
+   ! 1. Filling g%ind
+   ! The algorithm to get the folded index associated with a folded offset 
+   ! relies on the definition of the lexicographic order. 
+   ! Consideri appended upfront.ng α = (α₁,…,αₘ) with αᵢ ∈ { 1, ..., n },
+   ! the next index α' is such that there exists i that verifies
+   ! αⱼ = αⱼ' for all j < i, αᵢ' > αᵢ. Note that all the coordinates
+   ! αᵢ', ... , αₘ' need to be as small as the lexicographic order allows
+   ! for α' to immediately follow α.
+   ! Suppose j is the latest incremented coordinate: 
+   ! if αⱼ < n, then αⱼ' = αⱼ + 1
+   ! otherwise αⱼ = n, set αₖ' =  αⱼ₋₁ + 1 for all k ≥ j-1
+   ! if αⱼ₋₁ = n, set j := j-1  
+   ! otherwise, set j := m
+   ! 2. Filling g%count
+   ! The algorithm to count the number of equivalent unfolded indices
+   ! works as follows. A folded index can be written as 
+   ! α = (x₁, ..., x₁, ..., xₚ, ..., xₚ) such that x₁ < x₂ < ... < xₚ.
+   ! Denote kᵢ the number of coordinates equal to xᵢ.
+   ! The number of unfolded indices equivalent to α is 
+   ! c(α) = ⎛        m        ⎞
+   !        ⎝ k₁, k₂, ..., kₚ ⎠
+   ! k is an array such that k(ℓ,j) contains the number of coordinates
+   ! equal to ℓ for the folded index asociated with offset j. 
+   ! Suppose j is the latest incremented coordinate.
+   ! If αⱼ < n, then αⱼ' = αⱼ + 1, k(αⱼ) := k(αⱼ)-1, k(αⱼ') := k(αⱼ')+1.
+   ! In this case, c(α') = c(α)*(k(αⱼ)+1)/k(αⱼ')
+   ! otherwise, αⱼ = n: set αₖ' =  αⱼ₋₁ + 1 for all k ≥ j-1,
+   ! k(αⱼ₋₁) := k(αⱼ₋₁)-1, k(n) := 0, k(αⱼ₋₁') = m-(j-1)+1
+   ! In this case, we compute c(α') with the multinomial formula above
+   ! 3. Filling g%s
+   ! Suppose j is the latest incremented coordinate.
+   ! If αⱼ < n, then αⱼ' = αⱼ + 1: s' = s-1 if αⱼ = ns and s'=s otherwise
+   ! Otherwise, αⱼ = n: set αₖ' =  αⱼ₋₁ + 1 for all k ≥ j-1. Thus, 
+   ! s' = m if αⱼ₋₁ < ns; s'=s-1 if αⱼ₋₁ = ns; s'=s  otherwise
+   subroutine fill_folded_tensor(g, k, m, s, n, p)
+      type(folded_tensor), intent(inout) :: g
+      integer, contiguous, intent(inout) :: k(:,:)
+      integer, intent(in) :: m, s, n
+      type(pascal_triangle) :: p
+      integer :: j, lastinc
+      g%ind(1) = index(m, 1)
+      g%count(1) = 1
+      g%s(1) = m
+      k = 0
+      k(:,1) = 0
+      k(1,1) = m
+      j = 2
+      lastinc = m
+      do while (j <= size(g%ind))
+         g%ind(j) = index(g%ind(j-1)%coor)
+         k(:,j) = k(:,j-1)
+         if (g%ind(j-1)%coor(lastinc) == n) then
+            g%ind(j)%coor(lastinc-1:m) = g%ind(j-1)%coor(lastinc-1)+1
+            k(g%ind(j-1)%coor(lastinc-1),j) = k(g%ind(j-1)%coor(lastinc-1),j-1)-1
+            k(n,j) = 0
+            k(g%ind(j)%coor(lastinc-1),j) = m - (lastinc-1) + 1
+            g%count(j) = multinomial(k(:,j),m,p)
+            if (g%ind(j-1)%coor(lastinc-1) < s) then
+               g%s(j) = m 
+            elseif (g%ind(j-1)%coor(lastinc-1) == s) then
+               g%s(j) = g%s(j-1)-1
+            else
+               g%s(j) = g%s(j-1)
+            end if
+            if (g%ind(j)%coor(m) == n) then 
+               lastinc = lastinc-1
+            else
+               lastinc = m
+            end if
+         else
+            g%ind(j)%coor(lastinc) = g%ind(j-1)%coor(lastinc)+1
+            k(g%ind(j)%coor(lastinc),j) = k(g%ind(j)%coor(lastinc),j-1)+1
+            g%count(j) = g%count(j-1)*k(g%ind(j-1)%coor(lastinc),j)/k(g%ind(j)%coor(lastinc),j)
+            k(g%ind(j-1)%coor(lastinc),j) = k(g%ind(j-1)%coor(lastinc),j-1)-1
+            if (g%ind(j-1)%coor(lastinc) == s) then
+               g%s(j) = g%s(j-1)-1 
+            else
+               g%s(j) = g%s(j-1)
+            end if
+         end if
+         j = j+1
+      end do
+   end subroutine
+
+   ! Fills a tensor [gᵥᵐ] given the tensor [gᵥᵐ⁺¹] 
+   ! 1. Filling g%ind 
+   ! If (α₁,…,αₘ₊₁) is a folded index for [gᵥᵐ⁺¹], then (α₂,…,αₘ₊₁) is a folded
+   ! index of [gᵥᵐ]. The folded indices of [gᵥᵐ] are thus the tails of the first
+   ! ⎛ m+n-1 ⎞ folded indices of [gᵥᵐ⁺¹]
+   ! ⎝   m   ⎠
+   ! 2. Filling g%count
+   ! If α=(α₁,…,αₘ₊₁) is a folded index for [gᵥᵐ⁺¹], then α'=(α₂,…,αₘ₊₁) is a
+   ! folded index of [gᵥᵐ]. We thus have c(α') = c(α)*k(α₁)/(m+1) and 
+   ! perform k(α₁):=k(α₁)-1.
+   ! 3. Filling g%s
+   ! If α=(α₁,…,αₘ₊₁) is a folded index for [gᵥᵐ⁺¹], then α'=(α₂,…,αₘ₊₁) is a
+   ! folded index of [gᵥᵐ]. We thus have s(α') = s(α)-1 if α₁ ≤ s and 
+   ! s(α') = s(α) otherwise
+   subroutine fill_folded_tensor_backward(g_m, k, g_mp1, m, s)
+      type(folded_tensor), intent(inout) :: g_m ! tensor [gᵥᵐ]
+      ! k array from previous fill_folded_tensor_backward or 
+      ! fill_folded_tensor call. See definition in fill_folded_tensor
+      integer, contiguous, intent(inout) :: k(:,:)
+      type(folded_tensor), intent(in) :: g_mp1 ! tensor [gᵥᵐ⁺¹]
+      integer, intent(in) :: m, s ! s is the number of state variables
+      integer :: j
+      do j=1,size(g_m%ind)
+         g_m%ind(j)%coor = g_mp1%ind(j)%coor(2:m+1)
+         g_m%count(j) = g_mp1%count(j)*k(g_mp1%ind(j)%coor(1),j)/(m+1)
+         k(g_mp1%ind(j)%coor(1),j) = k(g_mp1%ind(j)%coor(1),j)-1
+         if (g_mp1%ind(j)%coor(1) <= s) then
+            g_m%s(j) = g_mp1%s(j)-1
+         else
+            g_m%s(j) = g_mp1%s(j)
+         end if
+      end do
+   end subroutine
+
+   ! Fills the ind, count and s members of a list of folded tensors g
+   ! with a number of state variables ns. More precisely, g contains the
+   ! information associated with the folded tensors [gᵥᵐ] with m=1,...,q
+   subroutine fill_list_folded_tensor(g, s, n, p)
+      type(folded_tensor), dimension(:), intent(inout) :: g
+      integer, intent(in) :: s, n ! number of state variables and size of v
+      type(pascal_triangle) :: p
+      integer :: q, m 
+      integer, dimension(:,:), allocatable :: k
+      q = size(g)
+      ! Case m = q
+      m = q
+      call allocate_folded_tensor(g(m), m, n, p)
+      allocate(k(n,size(g(m)%ind)))
+      call fill_folded_tensor(g(m), k, m, s, n, p) 
+      ! Case m < q
+      do m=q-1,1,-1
+         call allocate_folded_tensor(g(m), m, n, p)
+         call fill_folded_tensor_backward(g(m), k, g(m+1), m, s)
+      end do
+   end subroutine fill_list_folded_tensor
+end module tensors
+
+! After putting a comment on MATLAB-related lines
+! gfortran -o tensors tensors.f08 sort.f08 pascal.f08 partitions.f08
+! ./tensors
+! program test
+!    use pascal
+!    use tensors
+!    implicit none (type, external)
+!    type(pascal_triangle) :: p
+!    type(folded_tensor) :: g, h
+!    integer :: n, m, s, i, j
+!    integer, allocatable, dimension(:,:) :: k
+ 
+!    n = 3
+!    m = 3
+!    s = 2
+!    p = pascal_triangle(n+m-1)
+
+!    call allocate_folded_tensor(g, m, n, p)
+!    allocate(k(n, size(g%ind)))
+!    call fill_folded_tensor(g, k, c_null_ptr, m, s, n, p)
+
+!    ! List of folded indices, counts of equivalent unfolded indices
+!    ! and counts of state variables of [gᵥᵐ]
+!    ! Folded indices and offsets
+!    ! 0,0,0  1     1,1,1   7      2,2,2   10
+!    ! 0,0,1  2     1,1,2   8
+!    ! 0,0,2  3     1,2,2   9
+!    ! 0,1,1  4
+!    ! 0,1,2  5
+!    ! 0,2,2  6
+!    do i=1, size(k,2)
+!       print '(3i2)', (g%ind(i)%coor(j), j=1,m)
+!       print '(i2)', g%count(i)
+!       print '(i2)', g%s(i)
+!    end do
+
+!    call allocate_folded_tensor(h, m-1, n, p)
+!    call fill_folded_tensor_backward(h, k, g, c_null_ptr, m-1, s)
+!    ! List of folded indices, counts of equivalent unfolded indices
+!    ! and counts of state variables of [gᵥᵐ⁻¹]
+!    do i=1, size(h%ind)
+!       print '(2i2)', (h%ind(i)%coor(j), j=1,m-1)
+!       print '(i2)', h%count(i)
+!       print '(i2)', h%s(i)
+!    end do
+
+! end program test
\ No newline at end of file
diff --git a/mex/sources/local_state_space_iterations/local_state_space_iteration_k.f08 b/mex/sources/local_state_space_iterations/local_state_space_iteration_k.f08
index ce59cc47b5e7c82e8251f0247c43dd62a04b9e64..142aa34cad621ddbcde627926c6e333bb95366e5 100644
--- a/mex/sources/local_state_space_iterations/local_state_space_iteration_k.f08
+++ b/mex/sources/local_state_space_iterations/local_state_space_iteration_k.f08
@@ -26,7 +26,7 @@ module pparticle
    type tdata
       integer :: nm, nys, endo_nbr, nvar, order, nrestricted, nparticles 
       real(real64), allocatable :: yhat(:,:), e(:,:), ynext(:,:), ys_reordered(:), restrict_var_list(:)
-      type(pol), dimension(:), allocatable :: udr
+      type(tensor), dimension(:), allocatable :: udr
    end type tdata
 
    type(tdata) :: thread_data
@@ -37,7 +37,7 @@ contains
       type(c_ptr), intent(in), value :: arg
       integer, pointer :: im
       integer :: i, j, start, end, q, r, ind
-      type(horner), dimension(:), allocatable :: h
+      type(tensor), dimension(:), allocatable :: h
       real(real64), dimension(:), allocatable :: dyu
 
       ! Checking that the thread number got passed as argument
@@ -49,7 +49,7 @@ contains
       ! Allocating local arrays
       allocate(h(0:thread_data%order), dyu(thread_data%nvar)) 
       do i=0, thread_data%order
-         allocate(h(i)%c(thread_data%endo_nbr, thread_data%nvar**i))
+         allocate(h(i)%m(thread_data%endo_nbr, thread_data%nvar**i))
       end do
 
       ! Specifying bounds for the curent thread
@@ -69,7 +69,7 @@ contains
          call eval(h, dyu, thread_data%udr, thread_data%endo_nbr, thread_data%nvar, thread_data%order)
          do i=1,thread_data%nrestricted
             ind = int(thread_data%restrict_var_list(i))
-            thread_data%ynext(i,j) = h(0)%c(ind,1) + thread_data%ys_reordered(ind)
+            thread_data%ynext(i,j) = h(0)%m(ind,1) + thread_data%ys_reordered(ind)
          end do
       end do
 
@@ -100,7 +100,7 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
    type(c_ptr), dimension(*), intent(out) :: plhs
    integer(c_int), intent(in), value :: nlhs, nrhs
    type(c_ptr) :: M_mx, options_mx, dr_mx, yhat_mx, epsilon_mx, udr_mx, tmp
-   type(pol), dimension(:), allocatable :: udr
+   type(tensor), dimension(:), allocatable :: udr
    integer :: order, nstatic, npred, nboth, nfwrd, exo_nbr, endo_nbr, nparticles, nys, nvar, nrestricted, nm
    real(real64), dimension(:), pointer, contiguous :: order_var, ys, restrict_var_list
    real(real64), allocatable :: yhat(:,:), e(:,:), ynext(:,:), ys_reordered(:)
@@ -212,8 +212,8 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
       end if
       m = int(mxGetM(tmp))
       n = int(mxGetN(tmp))
-      allocate(udr(i)%g(m,n))
-      udr(i)%g(1:m,1:n) = reshape(mxGetPr(tmp), [m,n])
+      allocate(udr(i)%m(m,n))
+      udr(i)%m(1:m,1:n) = reshape(mxGetPr(tmp), [m,n])
    end do
 
    ! Initializing the global structure containing
diff --git a/tests/k_order_perturbation/fs2000_k_order_simul.mod b/tests/k_order_perturbation/fs2000_k_order_simul.mod
index a9133f0083de50d7d95f47cb118eb528faa2b23c..7c71fd4da5612739f505fbdd1fa9cc44651df9fe 100644
--- a/tests/k_order_perturbation/fs2000_k_order_simul.mod
+++ b/tests/k_order_perturbation/fs2000_k_order_simul.mod
@@ -67,22 +67,28 @@ steady;
 
 stoch_simul(order=2,nograph);
 
-ex=randn(10,M_.exo_nbr);
+ex=randn(5,M_.exo_nbr);
 
 Y2_matlab = simult_(M_,options_,oo_.dr.ys,oo_.dr,ex,options_.order);
 stoch_simul(order=2,k_order_solver);
-options_.k_order_solver=true;
 Y2_mex = simult_(M_,options_,oo_.dr.ys,oo_.dr,ex,options_.order); 
 
-if max(max(abs(Y2_matlab-Y2_mex)))>1e-8;
-   error('Matlab and mex simulation routines do not return same results')
+if max(abs(Y2_matlab(:)-Y2_mex(:)))>1e-8;
+   error('2nd order: Matlab and mex simulation routines do not return similar results')
 end
 
+Y2_mexiter = NaN(M_.endo_nbr,size(ex,1)+1);
+Y2_mexiter(:,1) = oo_.dr.ys;
+for it = 1:size(ex,1)
+    Y2_temp = simult_(M_,options_,Y2_mexiter(:,it),oo_.dr,ex(it,:),options_.order);
+    Y2_mexiter(:,it+1) = Y2_temp(:,2);    
+end
 
-stoch_simul(order=3,k_order_solver);
-
-ex=randn(5,M_.exo_nbr);
+if max((abs(Y2_mex(:) - Y2_mexiter(:))))>1e-8;
+   error('2nd order: sequential call does not return similar results')
+end
 
+stoch_simul(order=3, k_order_solver, nograph, irf=0);
 Y3_mex = simult_(M_,options_,oo_.dr.ys,oo_.dr,ex,options_.order); 
 
 Y3_mexiter = NaN(M_.endo_nbr,size(ex,1)+1);
@@ -92,6 +98,29 @@ for it = 1:size(ex,1)
     Y3_mexiter(:,it+1) = Y3_temp(:,2);    
 end
 
-if max(max(abs(Y3_mex - Y3_mexiter)))>1e-8;
-   error('mex simulation: sequential call does not return same results')
-end
\ No newline at end of file
+if max((abs(Y3_mex(:) - Y3_mexiter(:))))>1e-8;
+   error('3rd order: sequential call does not return similar results')
+end
+
+stoch_simul(order=2, k_order_solver, pruning, nograph, irf=0);
+
+options_.k_order_solver=false;
+Y2_matlab = simult_(M_,options_,oo_.dr.ys,oo_.dr,ex,options_.order); 
+options_.k_order_solver=true;
+Y2_mex = simult_(M_,options_,oo_.dr.ys,oo_.dr,ex,options_.order); 
+
+if max(abs(Y2_matlab(:)-Y2_mex(:)))>1e-8;
+   error('2nd order with pruning: Matlab and mex simulation routines do not return similar results')
+end
+
+stoch_simul(order=3, k_order_solver, pruning, nograph, irf=0);
+
+options_.k_order_solver=false;
+Y3_matlab = simult_(M_,options_,oo_.dr.ys,oo_.dr,ex,options_.order); 
+options_.k_order_solver=true;
+Y3_mex = simult_(M_,options_,oo_.dr.ys,oo_.dr,ex,options_.order); 
+
+if max(abs(Y3_matlab(:)-Y3_mex(:)))>1e-8;
+   error('3rd order with pruning: Matlab and mex simulation routines do not return similar results')
+end
+
diff --git a/tests/particle/first_spec_MCMC.mod b/tests/particle/first_spec_MCMC.mod
index 0bae888b5d7211e5efaf2aa2e6f426128436b4fe..9131f47c5d623a725ec81d4cc1dddcf8ad947c2c 100644
--- a/tests/particle/first_spec_MCMC.mod
+++ b/tests/particle/first_spec_MCMC.mod
@@ -15,7 +15,7 @@ stoch_simul(order=3,periods=200, irf=0);
 
 save('my_data_MCMC.mat','ca','b');
 
-options_.pruning=1;
+options_.pruning=true;
 options_.particle.pruning=true;
 options_.particle.number_of_particles=500;