diff --git a/mex/build/folded_to_unfolded_dr.am b/mex/build/folded_to_unfolded_dr.am
index c96e8a009..d0359723a 100644
--- a/mex/build/folded_to_unfolded_dr.am
+++ b/mex/build/folded_to_unfolded_dr.am
@@ -1,6 +1,6 @@
mex_PROGRAMS = folded_to_unfolded_dr
-folded_to_unfolded_dr_FCFLAGS = $(AM_FCFLAGS) -Warray-temporaries -I../libkordersim
+folded_to_unfolded_dr_FCFLAGS = $(AM_FCFLAGS) -I../libkordersim
nodist_folded_to_unfolded_dr_SOURCES = \
mexFunction.f08
diff --git a/mex/build/k_order_simul.am b/mex/build/k_order_simul.am
new file mode 100644
index 000000000..c7232168f
--- /dev/null
+++ b/mex/build/k_order_simul.am
@@ -0,0 +1,14 @@
+mex_PROGRAMS = k_order_simul
+
+k_order_simul_FCFLAGS = $(AM_FCFLAGS) -I../libkordersim
+
+nodist_k_order_simul_SOURCES = \
+ mexFunction.f08
+
+k_order_simul_LDADD = ../libkordersim/libkordersim.a
+
+BUILT_SOURCES = $(nodist_k_order_simul_SOURCES)
+CLEANFILES = $(nodist_k_order_simul_SOURCES)
+
+%.f08: $(top_srcdir)/../../sources/k_order_simul/%.f08
+ $(LN_S) -f $< $@
diff --git a/mex/build/local_state_space_iteration_fortran.am b/mex/build/local_state_space_iteration_fortran.am
index e3ee6eee7..1c7613394 100644
--- a/mex/build/local_state_space_iteration_fortran.am
+++ b/mex/build/local_state_space_iteration_fortran.am
@@ -1,6 +1,6 @@
mex_PROGRAMS = local_state_space_iteration_fortran
-local_state_space_iteration_fortran_FCFLAGS = $(AM_FCFLAGS) -Warray-temporaries -I../libkordersim
+local_state_space_iteration_fortran_FCFLAGS = $(AM_FCFLAGS) -I../libkordersim
nodist_local_state_space_iteration_fortran_SOURCES = \
mexFunction.f08
diff --git a/mex/build/matlab/Makefile.am b/mex/build/matlab/Makefile.am
index ae0ca758b..9609b1723 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 libkordersim folded_to_unfolded_dr local_state_space_iteration_fortran
+SUBDIRS += libdynare++ gensylv libkorder dynare_simul_ k_order_perturbation k_order_welfare local_state_space_iterations libkordersim folded_to_unfolded_dr local_state_space_iteration_fortran k_order_simul
endif
if ENABLE_MEX_MS_SBVAR
diff --git a/mex/build/matlab/configure.ac b/mex/build/matlab/configure.ac
index cfec90594..3476c4e6e 100644
--- a/mex/build/matlab/configure.ac
+++ b/mex/build/matlab/configure.ac
@@ -165,6 +165,7 @@ AC_CONFIG_FILES([Makefile
libkordersim/Makefile
folded_to_unfolded_dr/Makefile
local_state_space_iteration_fortran/Makefile
+ k_order_simul/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 000000000..99bfd2fb2
--- /dev/null
+++ b/mex/build/matlab/k_order_simul/Makefile.am
@@ -0,0 +1,2 @@
+include ../mex.am
+include ../../k_order_simul.am
\ No newline at end of file
diff --git a/mex/build/octave/Makefile.am b/mex/build/octave/Makefile.am
index 93295839b..53155961c 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 libkordersim folded_to_unfolded_dr local_state_space_iteration_fortran
+SUBDIRS += libdynare++ gensylv libkorder dynare_simul_ k_order_perturbation k_order_welfare local_state_space_iterations libkordersim folded_to_unfolded_dr local_state_space_iteration_fortran k_order_simul
endif
if ENABLE_MEX_MS_SBVAR
diff --git a/mex/build/octave/configure.ac b/mex/build/octave/configure.ac
index 7b4fa7b58..38ea69da4 100644
--- a/mex/build/octave/configure.ac
+++ b/mex/build/octave/configure.ac
@@ -158,6 +158,7 @@ AC_CONFIG_FILES([Makefile
libkordersim/Makefile
folded_to_unfolded_dr/Makefile
local_state_space_iteration_fortran/Makefile
+ k_order_simul/Makefile
perfect_foresight_problem/Makefile
num_procs/Makefile
block_trust_region/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 000000000..c214b9bc9
--- /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
\ No newline at end of file
diff --git a/mex/sources/Makefile.am b/mex/sources/Makefile.am
index 4a7ce297f..214c1382f 100644
--- a/mex/sources/Makefile.am
+++ b/mex/sources/Makefile.am
@@ -19,6 +19,7 @@ EXTRA_DIST = \
libkordersim \
folded_to_unfolded_dr \
local_state_space_iteration_fortran \
+ k_order_simul \
gensylv \
dynare_simul_ \
perfect_foresight_problem \
diff --git a/mex/sources/k_order_simul/mexFunction.f08 b/mex/sources/k_order_simul/mexFunction.f08
new file mode 100644
index 000000000..bf1f69ac6
--- /dev/null
+++ b/mex/sources/k_order_simul/mexFunction.f08
@@ -0,0 +1,180 @@
+! 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 .
+!
+! 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)
+! 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 struct
+ 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
+ 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
+ integer :: order, nstatic, npred, nboth, nfwrd, exo_nbr, endo_nbr, nys, nvar, nper
+ real(real64), dimension(:,:), allocatable :: shocks, sim
+ real(real64), dimension(:), allocatable :: ysteady, ystart, ysteady_pred, ystart_pred, dyu
+ type(pascal_triangle) :: p
+ type(uf_matching), dimension(:), allocatable :: matching
+ type(horner), dimension(:), allocatable :: h
+ integer :: i, t, d, m, n
+ character(kind=c_char, len=10) :: fieldname
+
+ order_mx = prhs(1)
+ nstatic_mx = prhs(2)
+ npred_mx = prhs(3)
+ nboth_mx = prhs(4)
+ nfwrd_mx = prhs(5)
+ nexog_mx = prhs(6)
+ ystart_mx = prhs(7)
+ shocks_mx = prhs(8)
+ ysteady_mx = prhs(9)
+ dr_mx = prhs(10)
+
+ ! 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")
+ end if
+ if (.not. (mxIsScalar(order_mx)) .and. mxIsNumeric(order_mx)) then
+ call mexErrMsgTxt("1st argument (order) should be a numeric scalar")
+ end if
+ if (.not. (mxIsScalar(nstatic_mx)) .and. mxIsNumeric(nstatic_mx)) then
+ call mexErrMsgTxt("2nd argument (nstat) should be a numeric scalar")
+ end if
+ if (.not. (mxIsScalar(npred_mx)) .and. mxIsNumeric(npred_mx)) then
+ call mexErrMsgTxt("3rd argument (npred) should be a numeric scalar")
+ end if
+ if (.not. (mxIsScalar(nboth_mx)) .and. mxIsNumeric(nboth_mx)) then
+ call mexErrMsgTxt("4th argument (nboth) should be a numeric scalar")
+ end if
+ if (.not. (mxIsScalar(nfwrd_mx)) .and. mxIsNumeric(nfwrd_mx)) then
+ call mexErrMsgTxt("5th argument (nforw) should be a numeric scalar")
+ end if
+ if (.not. (mxIsScalar(nexog_mx)) .and. mxIsNumeric(nexog_mx)) then
+ call mexErrMsgTxt("6th argument (nexog) should be a numeric scalar")
+ end if
+ if (.not. (mxIsDouble(ystart_mx) .and. (mxGetM(ystart_mx) == 1 .or. mxGetN(ystart_mx) == 1))) then
+ call mexErrMsgTxt("7th argument (ystart) should be a real vector")
+ end if
+ if (.not. (mxIsDouble(shocks_mx))) then
+ call mexErrMsgTxt("8th argument (shocks) should be a real matrix")
+ end if
+ if (.not. (mxIsDouble(ysteady_mx) .and. (mxGetM(ysteady_mx) == 1 .or. mxGetN(ysteady_mx) == 1))) then
+ call mexErrMsgTxt("9th argument (ysteady) should be a real vector")
+ end if
+ if (.not. mxIsStruct(dr_mx)) then
+ call mexErrMsgTxt("10th argument (dr) should be a struct")
+ end if
+
+ ! Converting inputs in Fortran format
+ order = int(mxGetScalar(order_mx))
+ nstatic = int(mxGetScalar(nstatic_mx))
+ npred = int(mxGetScalar(npred_mx))
+ nboth = int(mxGetScalar(nboth_mx))
+ nfwrd = int(mxGetScalar(nfwrd_mx))
+ exo_nbr = int(mxGetScalar(nexog_mx))
+ endo_nbr = nstatic+npred+nboth+nfwrd
+ nys = npred+nboth
+ nvar = nys+exo_nbr
+
+ if (endo_nbr /= int(mxGetM(ystart_mx))) then
+ call mexErrMsgTxt("ystart should have nstat+npred+nboth+nforw rows")
+ end if
+ allocate(ystart(endo_nbr))
+ ystart = mxGetPr(ystart_mx)
+
+ if (exo_nbr /= int(mxGetM(shocks_mx))) then
+ 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])
+
+ if (.not. (int(mxGetM(ysteady_mx)) == endo_nbr)) then
+ call mexErrMsgTxt("ysteady should have nstat+npred+nboth+nforw rows")
+ end if
+ allocate(ysteady(endo_nbr))
+ ysteady = mxGetPr(ysteady_mx)
+
+ 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))) then
+ call mexErrMsgTxt(trim(fieldname)//" is not allocated in dr")
+ 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)
+ 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)
+ udr(d)%g = fdr(d)%g(:,matching(d)%folded)
+ end do
+ 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)/))
+
+end subroutine mexFunction
\ No newline at end of file
diff --git a/mex/sources/local_state_space_iteration_fortran/mexFunction.f08 b/mex/sources/local_state_space_iteration_fortran/mexFunction.f08
index 81195ec83..a40a554a6 100644
--- a/mex/sources/local_state_space_iteration_fortran/mexFunction.f08
+++ b/mex/sources/local_state_space_iteration_fortran/mexFunction.f08
@@ -16,21 +16,14 @@
! along with Dynare. If not, see .
! 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,…)
-
+! yhat values of endogenous variables
+! epsilon values of the exgogenous shock
+! dr struct containing the folded tensors g_0, g_1, ...
+! M struct containing the model features
+! options struct containing the model options
+! udr struct containing the model unfolded tensors
! output:
-! res simulated results
+! ynext simulated next-period results
subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
use iso_fortran_env
@@ -45,14 +38,12 @@ 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, target :: fdr, udr
+ type(pol), dimension(:), allocatable, target :: udr
integer :: order, nstatic, npred, nboth, nfwrd, exo_nbr, endo_nbr, nparticles, nys, nvar, nrestricted
real(real64), dimension(:), allocatable :: order_var, ys, ys_reordered, restrict_var_list, dyu
real(real64), dimension(:,:), allocatable :: yhat, e, ynext, ynext_all
- type(pascal_triangle) :: p
- type(uf_matching), dimension(:), allocatable :: matching
type(horner), dimension(:), allocatable :: h
- integer :: i, d, j, m, n
+ integer :: i, j, m, n
character(kind=c_char, len=10) :: fieldname
yhat_mx = prhs(1)
@@ -63,9 +54,6 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
udr_mx = prhs(6)
! Checking the consistence and validity of input arguments
- ! if (nrhs /= 5 .or. nlhs /= 1) then
- ! call mexErrMsgTxt("Must have exactly 5 inputs and 1 output")
- ! end if
if (nrhs /= 6 .or. nlhs /= 1) then
call mexErrMsgTxt("Must have exactly 5 inputs and 1 output")
end if
@@ -142,7 +130,7 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
yhat = reshape(mxGetPr(yhat_mx), [nys, nparticles])
e = reshape(mxGetPr(epsilon_mx), [exo_nbr, nparticles])
- allocate(h(0:order), fdr(0:order), udr(0:order))
+ allocate(h(0:order), udr(0:order))
do i = 0, order
write (fieldname, '(a2, i1)') "g_", i
tmp = mxGetField(udr_mx, 1_mwIndex, trim(fieldname))
diff --git a/tests/k_order_perturbation/burnside_k_order.mod b/tests/k_order_perturbation/burnside_k_order.mod
index 000feeb0e..865ad0da0 100644
--- a/tests/k_order_perturbation/burnside_k_order.mod
+++ b/tests/k_order_perturbation/burnside_k_order.mod
@@ -123,3 +123,36 @@ 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'];
+ysteady = oo_.dr.ys(oo_.dr.order_var);
+dr = oo_.dr;
+
+vcov = M_.Sigma_e;
+seed = options_.DynareRandomStreams;
+
+tStart1 = tic; fortran_endo_simul = k_order_simul(order, nstat, npred, nboth, nfwrd, nexog, ystart, ex_, ysteady, dr); tElapsed1 = toc(tStart1);
+tStart2 = tic; dynare_endo_simul = dynare_simul_(order, nstat, npred, nboth, nfwrd, nexog, ystart,ex_,vcov,seed, ysteady, dr); tElapsed2 = toc(tStart2);
+
+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;
+
+if tElapsed1