diff --git a/matlab/utilities/tests b/matlab/utilities/tests
index 730078525361fb2896b04384ef4bc49a02a04c1c..bc4df1c1ec89a338424b5eb617fa568a76f14b45 160000
--- a/matlab/utilities/tests
+++ b/matlab/utilities/tests
@@ -1 +1 @@
-Subproject commit 730078525361fb2896b04384ef4bc49a02a04c1c
+Subproject commit bc4df1c1ec89a338424b5eb617fa568a76f14b45
diff --git a/mex/build/k_order_simul.am b/mex/build/k_order_simul.am
new file mode 100644
index 0000000000000000000000000000000000000000..b67a4532b824c3ebee1314275890e21df27db84e
--- /dev/null
+++ b/mex/build/k_order_simul.am
@@ -0,0 +1,28 @@
+mex_PROGRAMS = k_order_simul 
+
+nodist_k_order_simul_SOURCES = \
+	pascal.f08 \
+	sort.f08 \
+	partitions.f08 \
+	simulation.f08 \
+	mexFunction.f08 \
+	matlab_mex.F08 \
+	blas_lapack.F08
+
+BUILT_SOURCES = $(nodist_k_order_simul_SOURCES)
+CLEANFILES = $(nodist_k_order_simul_SOURCES)
+
+sort.mod: sort.o
+
+pascal.mod: pascal.o
+
+partitions.o: sort.mod pascal.mod
+partitions.mod : partitions.o
+
+simulation.o: partitions.mod 
+simulation.mod: simulation.o
+
+mexFunction.o: matlab_mex.mod simulation.mod partitions.mod pascal.mod
+
+%.f08: $(top_srcdir)/../../sources/k_order_simul/%.f08
+	$(LN_S) -f $< $@
diff --git a/mex/build/matlab/Makefile.am b/mex/build/matlab/Makefile.am
index 58f376cdc5e0614ec03d96f8e39c2909a302cd38..d33b580d716e16e1caf0e8d9ca7fb5247a6e7448 100644
--- a/mex/build/matlab/Makefile.am
+++ b/mex/build/matlab/Makefile.am
@@ -4,7 +4,7 @@ SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol perfect_foresight
 
 # libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_
 if ENABLE_MEX_DYNAREPLUSPLUS
-SUBDIRS += libdynare++ gensylv libkorder dynare_simul_ k_order_perturbation k_order_welfare local_state_space_iterations
+SUBDIRS += libdynare++ gensylv libkorder dynare_simul_ k_order_perturbation k_order_welfare k_order_simul local_state_space_iterations
 endif
 
 if ENABLE_MEX_MS_SBVAR
diff --git a/mex/build/matlab/configure.ac b/mex/build/matlab/configure.ac
index e838708a2b8b401495f4df15dbe9ab4599ec6b04..26f901eabcbf5937378893d0d758941c03294184 100644
--- a/mex/build/matlab/configure.ac
+++ b/mex/build/matlab/configure.ac
@@ -156,12 +156,13 @@ AC_CONFIG_FILES([Makefile
                  libkorder/Makefile
                  k_order_perturbation/Makefile
                  k_order_welfare/Makefile
+                 k_order_simul/Makefile
                  dynare_simul_/Makefile
                  kalman_steady_state/Makefile
                  ms_sbvar/Makefile
                  block_kalman_filter/Makefile
-	             sobol/Makefile
-		         local_state_space_iterations/Makefile
+                sobol/Makefile
+               local_state_space_iterations/Makefile
                  perfect_foresight_problem/Makefile
                  num_procs/Makefile
                  block_trust_region/Makefile
diff --git a/mex/build/matlab/k_order_simul/Makefile.am b/mex/build/matlab/k_order_simul/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..dea04644f572bdce84d0ef161d62a98ddd276ec6
--- /dev/null
+++ b/mex/build/matlab/k_order_simul/Makefile.am
@@ -0,0 +1,2 @@
+include ../mex.am
+include ../../k_order_simul.am
diff --git a/mex/build/octave/Makefile.am b/mex/build/octave/Makefile.am
index a669cf8393c3bd88866dc73c79f7f0dae6084df3..4daf7984ac63cc4e48bfa0a500f358cf95c7691c 100644
--- a/mex/build/octave/Makefile.am
+++ b/mex/build/octave/Makefile.am
@@ -4,7 +4,7 @@ SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol perfect_foresight
 
 # libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_
 if ENABLE_MEX_DYNAREPLUSPLUS
-SUBDIRS += libdynare++ gensylv libkorder dynare_simul_ k_order_perturbation k_order_welfare local_state_space_iterations
+SUBDIRS += libdynare++ gensylv libkorder dynare_simul_ k_order_perturbation k_order_welfare k_order_simul local_state_space_iterations
 endif
 
 if ENABLE_MEX_MS_SBVAR
diff --git a/mex/build/octave/configure.ac b/mex/build/octave/configure.ac
index 6a38305d059ec92d3a28ffa40cfdc03ff0dfbf90..cfdbcec108066a98a7e587335af6ea3e0923ce9a 100644
--- a/mex/build/octave/configure.ac
+++ b/mex/build/octave/configure.ac
@@ -138,6 +138,7 @@ AC_CONFIG_FILES([Makefile
                  libkorder/Makefile
                  k_order_perturbation/Makefile
                  k_order_welfare/Makefile
+                 k_order_simul/Makefile
                  dynare_simul_/Makefile
                  kalman_steady_state/Makefile
                  ms_sbvar/Makefile
diff --git a/mex/build/octave/k_order_simul/Makefile.am b/mex/build/octave/k_order_simul/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..186a19fe6e17677446cfd1e927b89a5998893c31
--- /dev/null
+++ b/mex/build/octave/k_order_simul/Makefile.am
@@ -0,0 +1,3 @@
+EXEEXT = .mex
+include ../mex.am
+include ../../k_order_simul.am
diff --git a/mex/sources/Makefile.am b/mex/sources/Makefile.am
index 1840bbf6a39724ae6782d6bd76c4858c6aad79e4..00edff535f9638d2773c62d7bb54c504ec8065f7 100644
--- a/mex/sources/Makefile.am
+++ b/mex/sources/Makefile.am
@@ -11,6 +11,7 @@ EXTRA_DIST = \
 	bytecode \
 	k_order_perturbation \
 	k_order_welfare \
+	k_order_simul \
 	kalman_steady_state \
 	ms-sbvar \
 	block_kalman_filter \
diff --git a/mex/sources/k_order_simul/mexFunction.f08 b/mex/sources/k_order_simul/mexFunction.f08
new file mode 100644
index 0000000000000000000000000000000000000000..7580b5a3775f0523c0b2fccf6c65fe3f933974ed
--- /dev/null
+++ b/mex/sources/k_order_simul/mexFunction.f08
@@ -0,0 +1,201 @@
+! Copyright © 2021 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/>.
+
+!  input:
+!       order    the order of approximation, needs order+1 derivatives
+!       nstat
+!       npred
+!       nboth
+!       nforw
+!       nexog
+!       ystart   starting value (full vector of endogenous)
+!       shocks   matrix of shocks (nexog x number of period)
+!       vcov     covariance matrix of shocks (nexog x nexog)
+!       seed     integer seed
+!       ysteady  full vector of decision rule's steady
+!       dr       structure containing matrices of derivatives (g_0, g_1,…)
+
+!  output:
+!       res      simulated results
+
+subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
+   use iso_fortran_env
+   use iso_c_binding
+   use matlab_mex
+   use partitions
+   use simulation
+   implicit none
+
+
+   type(c_ptr), dimension(*), intent(in), target :: prhs
+   type(c_ptr), dimension(*), intent(out) :: plhs
+   integer(c_int), intent(in), value :: nlhs, nrhs
+   real(real64), dimension(:), allocatable, target :: ystart, ysteady
+   real(real64), dimension(:,:), allocatable, target :: vcov, shocks, sim
+   type(pol), dimension(:), allocatable, target :: fdr
+   integer :: order, nstat, npred, nboth, nforw, nexog, seed, ny, nper, i, t 
+
+   ! Checking the consistence and validity of input arguments
+   if (nrhs /= 12 .or. nlhs /= 1) then
+      call mexErrMsgTxt("Must have exactly 12 inputs and 1 output")
+   end if
+   if (.not. (mxIsScalar(prhs(1))) .and. mxIsNumeric(prhs(1))) then
+      call mexErrMsgTxt("1st argument (order) should be a numeric scalar")
+   end if
+   if (.not. (mxIsScalar(prhs(2))) .and. mxIsNumeric(prhs(2))) then
+      call mexErrMsgTxt("2nd argument (nstat) should be a numeric scalar")
+   end if
+   if (.not. (mxIsScalar(prhs(3))) .and. mxIsNumeric(prhs(3))) then
+      call mexErrMsgTxt("3rd argument (npred) should be a numeric scalar")
+   end if
+   if (.not. (mxIsScalar(prhs(4))) .and. mxIsNumeric(prhs(4))) then
+      call mexErrMsgTxt("4th argument (nboth) should be a numeric scalar")
+   end if
+   if (.not. (mxIsScalar(prhs(5))) .and. mxIsNumeric(prhs(5))) then
+      call mexErrMsgTxt("5th argument (nforw) should be a numeric scalar")
+   end if
+   if (.not. (mxIsScalar(prhs(6))) .and. mxIsNumeric(prhs(6))) then
+      call mexErrMsgTxt("6th argument (nexog) should be a numeric scalar")
+   end if
+   if (.not. (mxIsDouble(prhs(7)) .and. (mxGetM(prhs(7)) == 1 .or. mxGetN(prhs(7)) == 1))) then
+      call mexErrMsgTxt("7th argument (ystart) should be a real vector")
+   end if
+   if (.not. (mxIsDouble(prhs(8)))) then
+      call mexErrMsgTxt("8th argument (shocks) should be a real matrix")
+   end if
+   if (.not. (mxIsDouble(prhs(9)))) then
+      call mexErrMsgTxt("9th argument (vcov) should be a real matrix")
+   end if
+   if (.not. (mxIsScalar(prhs(10))) .and. mxIsNumeric(prhs(10))) then
+      call mexErrMsgTxt("10th argument (seed) should be a numeric scalar")
+   end if
+   if (.not. (mxIsDouble(prhs(11)) .and. (mxGetM(prhs(11)) == 1 .or. mxGetN(prhs(11)) == 1))) then
+      call mexErrMsgTxt("11th argument (ysteady) should be a real vector")
+   end if
+   if (.not. mxIsStruct(prhs(12))) then
+      call mexErrMsgTxt("12th argument (dr) should be a struct")
+   end if
+
+   ! Converting inputs in Fortran format
+   order = int(mxGetScalar(prhs(1)))
+   nstat = int(mxGetScalar(prhs(2)))
+   npred = int(mxGetScalar(prhs(3)))
+   nboth = int(mxGetScalar(prhs(4)))
+   nforw = int(mxGetScalar(prhs(5)))
+   nexog = int(mxGetScalar(prhs(6)))
+   seed = int(mxGetScalar(prhs(10)))
+   ny = nstat + npred + nboth + nforw
+
+   associate (ystart_mat => mxGetPr(prhs(7)))
+      if (ny /= size(ystart_mat, 1)) then
+         call mexErrMsgTxt("ystart has wrong number of rows")
+      end if
+      allocate(ystart(size(ystart_mat)))
+      ystart = ystart_mat
+   end associate
+
+   associate (m => int(mxGetM(prhs(8))), n => int(mxGetN(prhs(8))), shocks_mat => mxGetPr(prhs(8)))
+      allocate(shocks(m, n))
+      shocks = reshape(shocks_mat, (/m, n/))
+      nper = n 
+   end associate
+
+   associate (m => int(mxGetM(prhs(9))), n => int(mxGetN(prhs(9))) , vcov_mat => mxGetPr(prhs(9)))
+      if (nexog /= m) then
+         call mexErrMsgTxt("vcov has a wrong number of rows")
+      end if
+      if (nexog /= n) then
+         call mexErrMsgTxt("vcov has a wrong number of cols")
+      end if
+      allocate(vcov(m,n))
+      vcov = reshape(vcov_mat, (/m,n/))
+   end associate
+
+   associate (ysteady_mat => mxGetPr(prhs(11)))
+      if (ny /= size(ysteady_mat, 1)) then
+         call mexErrMsgTxt("ysteady has wrong number of rows.\n")
+      end if
+      allocate(ysteady(size(ysteady_mat)))
+      ysteady = ysteady_mat
+   end associate
+
+   allocate(fdr(0:order)) 
+   do i = 0, order
+      block
+         type(c_ptr) :: tmp
+         character(len=10) :: str
+         integer :: m, n
+         write (str, '(a2, i1)') "g_", i
+         tmp = mxGetField(prhs(12), 1_mwIndex, trim(str))
+         if (.not. (c_associated(tmp) .and. mxIsDouble(tmp))) then
+            call mexErrMsgTxt(trim(str)//" is not allocated in dr")
+         end if
+         m = int(mxGetM(tmp))
+         n = int(mxGetN(tmp))
+         allocate(fdr(i)%g(m,n))
+         fdr(i)%g = reshape(mxGetPr(tmp), (/m, n/))
+      end block
+   end do
+
+   block
+      integer :: nys, nvar, d
+      real(real64), dimension(:), allocatable, target :: ystart_pred, ysteady_pred
+      type(pascal_triangle) :: p
+      type(uf_matching), dimension(:), allocatable :: matching 
+      real(real64), dimension (:,:), allocatable :: dyu, horner
+
+      allocate(sim(ny,nper))
+
+      ! Useful integers
+      nys = npred+nboth      ! number of pre-determined variables
+      nvar = nys+nexog       ! total number of "state" variables
+
+      if (order > 1) then
+         ! Compute the useful binomial coefficients from Pascal's triangle
+         p = pascal_triangle(nvar+order-1)
+         allocate(matching(2:order))
+         ! 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) 
+         end do
+      end if
+
+      allocate(dyu(nvar,1), ystart_pred(nys), ysteady_pred(nys))
+
+      ! Getting the predetermined part of the endogeneous variable vector 
+      ystart_pred = ystart(nstat+1:nstat+nys)
+      ysteady_pred = ysteady(nstat+1:nstat+nys) 
+      dyu(1:nys,1) = ystart_pred - ysteady_pred 
+      dyu(nys+1:,1) = shocks(:,1) 
+      call eval(horner, dyu, fdr, matching, ny, nvar, order)
+      sim(:,1) = horner(:,1) + ysteady
+     
+      ! Carrying out the simulation
+      do t=2,nper
+         dyu(1:nys,1) = horner(nstat+1:nstat+nys,1) 
+         dyu(nys+1:,1) = shocks(:,t)
+         call eval(horner, dyu, fdr, matching, ny, nvar, order)
+         sim(:,t) = horner(:,1) + ysteady
+      end do
+
+      plhs(1) = mxCreateDoubleMatrix(int(ny, mwSize), int(nper, mwSize), mxREAL)
+      mxGetPr(plhs(1)) = reshape(sim, (/size(sim)/))
+
+   end block
+
+end subroutine mexFunction
\ No newline at end of file
diff --git a/mex/sources/k_order_simul/partitions.f08 b/mex/sources/k_order_simul/partitions.f08
new file mode 100644
index 0000000000000000000000000000000000000000..861151b029a865444092fbe98adac331368fa2fd
--- /dev/null
+++ b/mex/sources/k_order_simul/partitions.f08
@@ -0,0 +1,276 @@
+! 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 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 partitions
+   use pascal
+   use sort
+   implicit none
+
+   ! index represents the aforementioned (α₁,…,αₘ) objects
+   type index
+      integer, dimension(:), allocatable :: ind 
+   end type index
+
+   interface index
+      module procedure :: init_index
+   end interface index
+
+   ! a dictionary that matches folded indices with folded offsets
+   type dict
+      integer :: pr ! pointer to the last added element in indices and offsets
+      type(index), dimension(:), allocatable :: indices ! list of folded indices
+      integer, dimension(:), allocatable :: offsets ! list of the associated offsets
+   end type dict
+
+   interface dict
+      module procedure :: init_dict
+   end interface dict 
+
+   interface operator(/=)
+      module procedure :: diff_indices
+   end interface operator(/=)
+
+   ! A type to contain the correspondence unfolded and folded offsets i.e.
+   ! folded(i) shall contain the folded offset corresponding to the unfolded offset i
+   type uf_matching
+      type(integer), dimension(:), allocatable :: folded
+   end type
+
+   ! ! eq_class defines the equivalent class of a given index
+   ! ! i.e. (1,2,1), (2,1,1) and (1,1,2) are equivalent indices
+   ! ! with representative (1,1,2)
+   ! type eq_class
+   !    type(index) :: r ! unique representative of the equivalent class
+   !    type(index), dimension(:), allocatable :: c ! class
+   ! end type eq_class
+
+   ! ! partition is meant to contain all the indexes associated with an
+   ! ! integer d
+   ! type partition
+   !    integer :: d ! length of the partition
+   !    integer :: n ! number of values αᵢ may take
+   !    type(index), dimension(:), allocatable :: p ! partition of indexes
+   ! end type partition
+   
+   ! interface partition
+   !    module procedure :: init_partition
+   ! end interface partition
+
+contains
+
+   ! Constructor for the index type
+   type(index) function init_index(d, ind)
+      integer, intent(in) :: d
+      integer, dimension(d), intent(in) :: ind
+      allocate(init_index%ind(d))
+      init_index%ind = ind
+   end function init_index
+
+   ! Comparison for the index type. Returns true if the two indices are different
+   type(logical) function diff_indices(i1,i2)
+      type(index), intent(in) :: i1, i2 
+      if (size(i1%ind) /= size(i2%ind) .or. any(i1%ind /= i2%ind)) then
+         diff_indices = .true.
+      else
+         diff_indices = .false.
+      end if
+   end function diff_indices
+
+   ! Constructor of the dict type
+   type(dict) function init_dict(n, d, p)
+      integer, intent(in) :: n, d
+      type(pascal_triangle), intent(in) :: p
+      integer :: size
+      size = get(d, n+d-1, p)
+      allocate(init_dict%indices(size), init_dict%offsets(size))
+      init_dict%pr = 0
+   end function init_dict
+
+   ! Count the number of coordinates similar to the the first one for a given index
+   type(integer) function get_prefix_length(idx, d)
+      integer, intent(in) :: d
+      type(index), intent(in) :: idx
+      integer :: i
+      i = 1
+      if (d>1) then
+         do while ((i < d) .and. (idx%ind(i+1) == idx%ind(1)))
+            i = i+1
+         end do
+      end if
+      get_prefix_length = i
+   end function get_prefix_length
+
+   ! Gets the folded index associated with an unfolded index
+   type(index) function u_index_to_f_index(idx, d)
+      type(index), intent(in) :: idx
+      integer, intent(in) :: d
+      u_index_to_f_index = index(d, idx%ind)
+      call sort_int(u_index_to_f_index%ind, d)
+   end function u_index_to_f_index 
+
+   ! Converts the offset of an unfolded tensor to the associated unfolded tensor index
+   ! Note that the index (α₁,…,αₘ) is such that αᵢ ∈ { 0, ..., n-1 }
+   ! and the offset is such that j ∈ {1, ..., nᵈ}
+   type(index) function u_offset_to_u_index(j, n, d)
+      integer, intent(in) :: j, n, d ! offset, number of variables and dimensions respectively
+      integer :: i, tmp, r
+      allocate(u_offset_to_u_index%ind(d))
+      tmp = j-1 ! We substract 1 as j ∈ {1, ..., n} so that tmp ∈ {0, ..., n-1} and our modular operations work
+      do i=d,1,-1
+         r = mod(tmp, n)
+         u_offset_to_u_index%ind(i) = r
+         tmp = (tmp-r)/n
+      end do
+   end function u_offset_to_u_index
+
+   ! Converts a folded tensor index to the associated folded tensor offset
+   ! See the explanation in dynare++/tl/cc/tensor.cc for the function FTensor::getOffsetRecurse
+   ! Note that the index (α₁,…,αₘ) is such that αᵢ ∈ { 0, ..., n-1 }
+   ! and the offset is such that j ∈ {1, ...,  ⎛n+d-1⎞ }
+   !                                           ⎝  d  ⎠
+   recursive function f_index_to_f_offset(idx, n, d, p) result(j)
+      type(index), intent(in) :: idx   ! folded index
+      integer, intent(in) :: n, d      ! number of variables and dimensions
+      type(pascal_triangle) :: p       ! Pascal's triangle containing the relevant binomial coefficients
+      integer :: j, prefix
+      type(index) :: tmp
+      if (d == 0) then
+         j = 1
+      else
+         prefix = get_prefix_length(idx,d)
+         tmp = index(d-prefix, idx%ind(prefix+1:) - idx%ind(1))
+         j = get(d, n+d-1, p) - get(d, n-idx%ind(1)+d-1, p) + f_index_to_f_offset(tmp, n-idx%ind(1), d-prefix, p) 
+      end if
+   end function f_index_to_f_offset
+ 
+   ! Function that searches a value in an array of a given length 
+   type(integer) function find(a, v, l)
+      integer, intent(in) :: l ! length of the array
+      type(index), dimension(l), intent(in) :: a ! array of indices
+      type(index) :: v ! element to be found
+      integer :: i
+      if (l == 0) then
+         find = 0
+      else
+         i = 1
+         do while (i <= l .and. a(i) /= v)
+            i = i+1
+         end do
+         if (i == l+1) then
+            find = 0
+         else
+            find = i
+         end if
+      end if
+   end function find
+
+   ! Fills the folded offset array:
+   ! folded(i) shall contain the folded offset corresponding to the unfolded offset i
+   ! For each unfolded tensor offset
+   ! 	(a) compute the associated unfolded index (u_offset_to_u_index)
+   ! 	(b) compute the associated folded index (u_index_to_f_index)
+   ! 	(c) has the folded offset already been computed ?
+   ! 		(i)   If yes, get the corresponding offset
+   ! 		(ii)  If no, compute it (f_index_to_f_offset) and store it for reuse and as a result
+   subroutine fill_folded_indices(folded, n, d, p) 
+      integer, intent(in) :: n, d
+      integer, dimension(n**d), intent(inout) :: folded 
+      type(pascal_triangle), intent(in) :: p
+      type(dict) :: c
+      type(index) :: tmp
+      integer :: j, found
+      c = dict(n, d, p)
+      do j=1,n**d
+         tmp = u_offset_to_u_index(j,n,d)
+         tmp = u_index_to_f_index(tmp, d)
+         found = find(c%indices, tmp, c%pr)
+         if (found == 0) then
+           c%pr = c%pr+1
+           c%indices(c%pr) = tmp
+           c%offsets(c%pr) = f_index_to_f_offset(tmp,n,d,p)
+           folded(j) = c%offsets(c%pr)  
+         else
+           folded(j) = c%offsets(found) 
+         end if
+      end do
+   end subroutine fill_folded_indices
+
+   end module partitions
+
+! gfortran -o partitions partitions.f08 pascal.f08 sort.f08
+! ./partitions
+! program test
+!    use partitions
+!    use pascal
+!    implicit none
+!    type(index) :: uidx, fidx, i1, i2
+!    integer, dimension(:), allocatable :: folded
+!    integer :: i, uj, n, d
+!    type(pascal_triangle) :: p
+!    ! Unfolded indices and offsets
+!    ! 0,0,0  1    1,0,0  10   2,0,0  19
+!    ! 0,0,1  2    1,0,1  11   2,0,1  20
+!    ! 0,0,2  3    1,0,2  12   2,0,2  21
+!    ! 0,1,0  4    1,1,0  13   2,1,0  22
+!    ! 0,1,1  5    1,1,1  14   2,1,1  23
+!    ! 0,1,2  6    1,1,2  15   2,1,2  24
+!    ! 0,2,0  7    1,2,0  16   2,2,0  25
+!    ! 0,2,1  8    1,2,1  17   2,2,1  26
+!    ! 0,2,2  9    1,2,2  18   2,2,2  27  
+
+!    ! 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
+ 
+!    n = 3
+!    d = 3
+!    uj = 8
+!    p = pascal_triangle(n+d-1)   
+
+!    ! u_offset_to_u_index
+!    uidx = u_offset_to_u_index(uj,n,d)
+!    print '(3i2)', (uidx%ind(i), i=1,d) ! should display 0 2 1
+
+!    ! f_index_to_f_offset
+!    fidx = u_index_to_f_index(uidx, d)
+!    print '(i2)', f_index_to_f_offset(fidx, n, d, p) ! should display 5
+
+!    ! /=
+!    i1 = index(5, (/1,2,3,4,5/))
+!    i2 = index(5, (/1,2,3,4,6/))
+!    if (i1 /= i2) then
+!       print *, "Same!"
+!    else
+!       print *, "Different!"
+!    end if
+
+!    ! fill_folded_indices
+!    allocate(folded(n**d))
+!    call fill_folded_indices(folded,n,d)
+!    print *, "Matching offsets unfolded -> folded"
+!    print '(1000i4)', (i, i=1,n**d)
+!    print '(1000i4)', (folded(i), i=1,n**d)
+
+! end program test
\ No newline at end of file
diff --git a/mex/sources/k_order_simul/pascal.f08 b/mex/sources/k_order_simul/pascal.f08
new file mode 100644
index 0000000000000000000000000000000000000000..bec24515ec0b50d12973886f359d1d8577cfafed
--- /dev/null
+++ b/mex/sources/k_order_simul/pascal.f08
@@ -0,0 +1,89 @@
+! Provides a subroutine to build Pascal's triangle
+!
+! Copyright © 2021 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 pascal
+   implicit none
+
+   type line
+      integer, dimension(:), allocatable :: coeffs 
+   end type line
+
+   type pascal_triangle
+      integer :: d
+      type(line), dimension(:), allocatable :: lines     
+   end type pascal_triangle
+
+   interface pascal_triangle
+      module procedure :: init_pascal_triangle
+   end interface pascal_triangle
+
+contains
+
+   ! Fills a pascal_triangle structure with the associated coefficients up to a given order
+   type(pascal_triangle) function init_pascal_triangle(d)
+      integer, intent(in) :: d
+      integer :: i, j
+      init_pascal_triangle%d = d
+      allocate(init_pascal_triangle%lines(d))
+      ! Initializing the first line
+      allocate(init_pascal_triangle%lines(1)%coeffs(2))
+      init_pascal_triangle%lines(1)%coeffs = 1
+      ! Iterating Pascal's triangle formula
+      if (d > 1) then
+         do i=2,d
+            allocate(init_pascal_triangle%lines(i)%coeffs(i+1))
+            init_pascal_triangle%lines(i)%coeffs(1) = 1
+            init_pascal_triangle%lines(i)%coeffs(i+1) = 1
+            do j=2,i
+               init_pascal_triangle%lines(i)%coeffs(j) = init_pascal_triangle%lines(i-1)%coeffs(j-1) &
+                                                      + init_pascal_triangle%lines(i-1)%coeffs(j)
+            end do
+         end do
+      end if
+   end function init_pascal_triangle
+
+   ! Returns ⎛n⎞ stored in pascal_triangle p
+   !         ⎝k⎠
+ 
+   type(integer) function get(k,n,p)
+      integer, intent(in) :: k, n
+      type(pascal_triangle), intent(in) :: p
+      get = p%lines(n)%coeffs(k+1)
+   end function get
+
+end module pascal
+
+! gfortran -o pascal pascal.f08
+! ./pascal
+! program test
+!    use pascal
+!    type(pascal_triangle) :: p
+!    integer :: d
+!    read *, d
+!    p = pascal_triangle(d)
+!    do n=1,d 
+!       do k=0,n
+!          if (k < n) then
+!             write (*,'(i2," ")', advance="no") get(k,n,p)
+!          else
+!             write (*,'(i2," ")') get(k,n,p)
+!          end if
+!       end do
+!    end do
+! end program test
\ No newline at end of file
diff --git a/mex/sources/k_order_simul/simulation.f08 b/mex/sources/k_order_simul/simulation.f08
new file mode 100644
index 0000000000000000000000000000000000000000..001703b589c27fcfa455eb956ef90e407838bd6d
--- /dev/null
+++ b/mex/sources/k_order_simul/simulation.f08
@@ -0,0 +1,74 @@
+! Necessary routines and functions to carry out simulations
+!
+! A first step is to get the associated 
+! Copyright © 2021 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 simulation
+   use iso_fortran_env
+   use partitions
+   implicit none
+
+   type :: pol
+      real(real64), dimension(:,:), allocatable :: g
+   end type pol
+
+contains
+
+   ! Horner evaluation of a tensor polynomial at a given vector
+   subroutine eval(horner, dyu, fdr, matching, ny, nvar, order)
+      real(real64), dimension(:,:), allocatable, intent(inout) :: horner
+      real(real64), dimension(nvar,1), intent(in) :: dyu
+      type(pol), intent(in) :: fdr(0:order)
+      type(uf_matching), intent(in) :: matching(2:order)
+      integer, intent(in) :: ny, nvar, order
+      integer :: d
+      if (order == 1) then
+         horner = fdr(1)%g
+      else
+         horner = contract(fdr(order)%g(:, matching(order)%folded), dyu, ny, nvar, order)
+      end if
+      do d=order-1,1,-1
+         if (d > 1) then
+            horner = horner + fdr(d)%g(:,matching(d)%folded)
+            horner = contract(horner, dyu, ny, nvar, d)
+         else
+            horner = horner + fdr(1)%g
+         end if
+      end do
+      horner = matmul(horner, dyu) + fdr(0)%g
+   end subroutine eval
+
+   ! Contracts the unfolded tensor t with respect to the vector x and stores the result in c
+   function contract(t, x, nrows, nvar, d)
+      integer, intent(in) :: nrows, nvar, d
+      real(real64), dimension(nrows, nvar**d), intent(in) :: t
+      real(real64), dimension(nvar,1), intent(in) :: x
+      real(real64), dimension(nrows, nvar**(d-1)) :: contract
+      integer :: i
+      contract = 0.
+      do i=1,nvar**(d-1)
+         contract(:,i) = reshape(matmul(t(:,(i-1)*nvar+1:i*nvar), x), (/nrows/))
+      end do
+   end function contract
+
+end module simulation
+
+! program test
+!    implicit none
+   
+! end program test
\ No newline at end of file
diff --git a/mex/sources/k_order_simul/sort.f08 b/mex/sources/k_order_simul/sort.f08
new file mode 100644
index 0000000000000000000000000000000000000000..9d2b02c736fe76b90344867b41c0e9b77dd309dc
--- /dev/null
+++ b/mex/sources/k_order_simul/sort.f08
@@ -0,0 +1,39 @@
+! Provides a subroutine to sort integer arrays in ascending order
+!
+! Copyright © 2021 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 sort
+   implicit none
+
+contains
+   subroutine sort_int(l, d)
+      integer, intent(in) :: d
+      integer, dimension(d), intent(inout) :: l
+      integer :: i, j, x
+      do i=2,d
+         x = l(i)
+         j = i
+         do while (j > 1 .and. l(j-1) > x)
+            l(j) = l(j-1)
+            j = j-1
+         end do
+         l(j) = x
+      end do
+   end subroutine sort_int
+
+end module sort
\ No newline at end of file
diff --git a/tests/k_order_perturbation/burnside_k_order.mod b/tests/k_order_perturbation/burnside_k_order.mod
index 000feeb0ebc00fce4fef103632dab97f3bf532a4..b7a642f39d14b83694fc203576ce2793d754d178 100644
--- a/tests/k_order_perturbation/burnside_k_order.mod
+++ b/tests/k_order_perturbation/burnside_k_order.mod
@@ -53,7 +53,7 @@
 
  */
 
-@#define k = 9
+@#define k = 3
 
 var y x;
 varexo e;
@@ -123,3 +123,23 @@ for T = 1:size(oo_.endo_simul,2)
   end
   xlag = oo_.endo_simul(2,T);
 end
+
+% Verify that the simulated time series is correct with the Fortran routine k_order_simul
+
+order = options_.order;
+nstat = M_.nstatic;
+npred = M_.npred;
+nboth = M_.nboth;
+nfwrd = M_.nfwrd;
+nexog = M_.exo_nbr;
+ystart = oo_.dr.ys(oo_.dr.order_var,1);
+ex_ = [zeros(M_.maximum_lag,M_.exo_nbr), oo_.exo_simul'];
+vcov = M_.Sigma_e;
+seed = options_.DynareRandomStreams;
+ysteady = oo_.dr.ys(oo_.dr.order_var);
+dr = oo_.dr;
+fortran_endo_simul = k_order_simul(order, nstat, npred, nboth, nfwrd, nexog, ystart, ex_, vcov, seed, ysteady, dr);
+
+if (max(abs(oo_.endo_simul-fortran_endo_simul(oo_.dr.order_var,2:end))) > 1e-10)
+    error('Error in k_order_simul: inaccurate simulation');
+end;
\ No newline at end of file