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