stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Adds pruning in k_order_simul 

See issue #1643 about beyond-third-order pruning
kalman-mex
Normann Rion 2023-07-05 14:34:07 +02:00
parent c02e550582
commit 710589eb5b
17 changed files with 836 additions and 149 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
matlab/k_order_pert.m
View File

@ -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)

9
matlab/simult_.m
View File

@ -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

6
mex/build/libkordersim.am
View File

@ -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

14
mex/sources/folded_to_unfolded_dr/mexFunction.f08
View File

@ -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

20
mex/sources/k_order_mean/mexFunction.f08
View File

@ -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

42
mex/sources/k_order_perturbation/k_order_perturbation.cc
View File

@ -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)

93
mex/sources/k_order_simul/mexFunction.f08
View File

@ -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

7
mex/sources/k_order_welfare/k_order_welfare.cc
View File

@ -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();

24
mex/sources/libkorder/kord/approximation.cc
View File

@ -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

6
mex/sources/libkorder/kord/approximation.hh
View File

@ -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

53
mex/sources/libkorder/kord/decision_rule.hh
View File

@ -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)

119
mex/sources/libkordersim/partitions.f08
View File

@ -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

219
mex/sources/libkordersim/simulation.f08
View File

@ -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

304
mex/sources/libkordersim/tensors.f08 Normal file
View File

@ -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

14
mex/sources/local_state_space_iterations/local_state_space_iteration_k.f08
View File

@ -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

49
tests/k_order_perturbation/fs2000_k_order_simul.mod
View File

@ -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
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

2
tests/particle/first_spec_MCMC.mod
View File

@ -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;