diff --git a/mex/sources/local_state_space_iterations/local_state_space_iteration_k.cc b/mex/sources/local_state_space_iterations/local_state_space_iteration_k.cc
index 43e9a4779f109951626f3fca6bf91ad6a80a93d3..d370cb36ed862a23b33907cc7168ae8b3683f715 100644
--- a/mex/sources/local_state_space_iterations/local_state_space_iteration_k.cc
+++ b/mex/sources/local_state_space_iterations/local_state_space_iteration_k.cc
@@ -1,5 +1,5 @@
/*
- * Copyright © 2019 Dynare Team
+ * Copyright © 2019-2021 Dynare Team
*
* This file is part of Dynare.
*
@@ -41,15 +41,16 @@ struct ParticleWorker : public sthread::detach_thread
const ConstGeneralMatrix &yhat, ε
const Vector &ys_reordered;
const UnfoldDecisionRule &dr;
+ const ConstVector &restrict_var_list;
GeneralMatrix &ynext;
ParticleWorker(int npred_both_arg, int exo_nbr_arg, std::pair<size_t, size_t> particle_range_arg,
const ConstGeneralMatrix &yhat_arg, const ConstGeneralMatrix &epsilon_arg,
const Vector &ys_reordered_arg, const UnfoldDecisionRule &dr_arg,
- GeneralMatrix &ynext_arg)
+ const ConstVector &restrict_var_list_arg, GeneralMatrix &ynext_arg)
: npred_both{npred_both_arg}, exo_nbr{exo_nbr_arg}, particle_range{std::move(particle_range_arg)},
yhat{yhat_arg}, epsilon{epsilon_arg}, ys_reordered{ys_reordered_arg}, dr{dr_arg},
- ynext{ynext_arg}
+ restrict_var_list{restrict_var_list_arg}, ynext{ynext_arg}
{
}
void
@@ -58,16 +59,22 @@ struct ParticleWorker : public sthread::detach_thread
Vector dyu(npred_both+exo_nbr);
Vector dy(dyu, 0, npred_both);
Vector u(dyu, npred_both, exo_nbr);
+ Vector ynext_col_allvars(ys_reordered.length());
for (size_t i = particle_range.first; i < particle_range.second; i++)
{
dy = yhat.getCol(i);
u = epsilon.getCol(i);
- Vector ynext_col{ynext.getCol(i)};
- dr.eval(DecisionRule::emethod::horner, ynext_col, dyu);
+ dr.eval(DecisionRule::emethod::horner, ynext_col_allvars, dyu);
+
+ ynext_col_allvars.add(1.0, ys_reordered);
- ynext_col.add(1.0, ys_reordered);
+ /* Select only the variables in restrict_var_list, and copy back to the
+ result matrix */
+ Vector ynext_col{ynext.getCol(i)};
+ for (int j = 0; j < restrict_var_list.length(); j++)
+ ynext_col[j] = ynext_col_allvars[restrict_var_list[j]-1];
}
}
};
@@ -107,6 +114,9 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
const mxArray *ys_mx = mxGetField(dr_mx, 0, "ys");
if (!ys_mx || !mxIsDouble(ys_mx) || mxGetNumberOfElements(ys_mx) != static_cast<size_t>(endo_nbr))
mexErrMsgTxt("Field dr.ys should be a double precision vector with endo_nbr elements");
+ const mxArray *restrict_var_list_mx = mxGetField(dr_mx, 0, "restrict_var_list");
+ if (!(restrict_var_list_mx && mxIsDouble(restrict_var_list_mx)))
+ mexErrMsgTxt("Field dr.restrict_var_list should be a double precision vector");
size_t nparticles = mxGetN(yhat_mx);
if (mxGetN(epsilon_mx) != nparticles)
@@ -125,6 +135,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
ConstGeneralMatrix epsilon{epsilon_mx};
ConstVector ys{ys_mx};
const double *order_var = mxGetPr(order_var_mx);
+ ConstVector restrict_var_list{restrict_var_list_mx};
try
{
@@ -163,7 +174,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
exo_nbr, ys_reordered);
// Create the result matrix
- plhs[0] = mxCreateDoubleMatrix(endo_nbr, nparticles, mxREAL);
+ plhs[0] = mxCreateDoubleMatrix(restrict_var_list.length(), nparticles, mxREAL);
GeneralMatrix ynext{plhs[0]};
// Run the real job in parallel
@@ -174,6 +185,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
for (size_t i = 0; i < nparticles; i += part_by_thread)
group.insert(std::make_unique<ParticleWorker>(npred+nboth, exo_nbr,
std::make_pair(i, std::min(i+part_by_thread, nparticles)),
- yhat, epsilon, ys_reordered, dr, ynext));
+ yhat, epsilon, ys_reordered, dr, restrict_var_list,
+ ynext));
group.run();
}
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 14e6c86f24a2798b7347a5f661d1ea8d059a4f62..0040c4b02f08756ee9fc46db7092bd1134f5093e 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -825,6 +825,10 @@ discretionary_policy/dennis_1_estim.o.trs: discretionary_policy/dennis_1.o.trs
discretionary_policy/Gali_2015_chapter_3_nonlinear.m.trs: discretionary_policy/Gali_2015_chapter_3.m.trs
discretionary_policy/Gali_2015_chapter_3_nonlinear.o.trs: discretionary_policy/Gali_2015_chapter_3.o.trs
+particle/local_state_space_iteration_k_test.m.trs: particle/first_spec.m.trs
+particle/local_state_space_iteration_k_test.o.trs: particle/first_spec.o.trs
+
+
observation_trends_and_prefiltering/MCMC: m/observation_trends_and_prefiltering/MCMC o/observation_trends_and_prefiltering/MCMC
m/observation_trends_and_prefiltering/MCMC: $(patsubst %.mod, %.m.trs, $(filter observation_trends_and_prefiltering/MCMC/%.mod, $(MODFILES)))
o/observation_trends_and_prefiltering/MCMC: $(patsubst %.mod, %.o.trs, $(filter observation_trends_and_prefiltering/MCMC/%.mod, $(MODFILES)))
diff --git a/tests/particle/local_state_space_iteration_k_test.mod b/tests/particle/local_state_space_iteration_k_test.mod
index 5c805344f70b899d4ad95cdecad6300627a33a74..bf09dc5017b82e9cc19a44a164341e5b0e3fd337 100644
--- a/tests/particle/local_state_space_iteration_k_test.mod
+++ b/tests/particle/local_state_space_iteration_k_test.mod
@@ -1,51 +1,29 @@
/*
Tests that local_state_space_iteration_2 and local_state_space_iteration_k
(for k=2) return the same results.
-*/
-
-var y, k, a, h, b, c;
-varexo e, u;
-parameters beta, rho, alpha, delta, theta, psi, tau;
+ This file must be run after first_spec.mod (both are based on the same model).
+*/
-alpha = 0.36;
-rho = 0.95;
-tau = 0.025;
-beta = 0.99;
-delta = 0.025;
-psi = 0;
-theta = 2.95;
+@#include "first_spec_common.inc"
-phi = 0.1;
+varobs q ca;
-model;
-c*theta*h^(1+psi)=(1-alpha)*y;
-k = beta*(((exp(b)*c)/(exp(b(+1))*c(+1)))
- *(exp(b(+1))*alpha*y(+1)+(1-delta)*k));
-y = exp(a)*(k(-1)^alpha)*(h^(1-alpha));
-k = exp(b)*(y-c)+(1-delta)*k(-1);
-a = rho*a(-1)+tau*b(-1) + e;
-b = tau*a(-1)+rho*b(-1) + u;
+shocks;
+var eeps = 0.04^2;
+var nnu = 0.03^2;
+var q = 0.01^2;
+var ca = 0.01^2;
end;
-initval;
-y = 1.08068253095672;
-c = 0.80359242014163;
-h = 0.29175631001732;
-k = 11.08360443260358;
-a = 0;
-b = 0;
-e = 0;
-u = 0;
-end;
+// Initialize various structures
+estimation(datafile='my_data.mat',order=2,mode_compute=0,mh_replic=0,filter_algorithm=sis,nonlinear_filter_initialization=2
+ ,cova_compute=0 %tell program that no covariance matrix was computed
+);
-shocks;
-var e; stderr 0.009;
-var u; stderr 0.009;
-var e, u = phi*0.009*0.009;
-end;
+stoch_simul(order=2, periods=200, irf=0, k_order_solver);
-stoch_simul(order=2, irf=0, k_order_solver);
+// Really perform the test
nparticles = 100;
@@ -57,7 +35,16 @@ epsilon = chol(M_.Sigma_e)*randn(M_.exo_nbr, nparticles);
dr = oo_.dr;
-tStart1 = tic; for i=1:10000, ynext1 = local_state_space_iteration_2(yhat, epsilon, dr.ghx, dr.ghu, dr.ys(dr.order_var)+0.5*dr.ghs2, dr.ghxx, dr.ghuu, dr.ghxu, 1); end, tElapsed1 = toc(tStart1);
+// “rf” stands for “Reduced Form”
+rf_ghx = dr.ghx(dr.restrict_var_list, :);
+rf_ghu = dr.ghu(dr.restrict_var_list, :);
+rf_constant = dr.ys(dr.order_var)+0.5*dr.ghs2;
+rf_constant = rf_constant(dr.restrict_var_list, :);
+rf_ghxx = dr.ghxx(dr.restrict_var_list, :);
+rf_ghuu = dr.ghuu(dr.restrict_var_list, :);
+rf_ghxu = dr.ghxu(dr.restrict_var_list, :);
+
+tStart1 = tic; for i=1:10000, ynext1 = local_state_space_iteration_2(yhat, epsilon, rf_ghx, rf_ghu, rf_constant, rf_ghxx, rf_ghuu, rf_ghxu, 1); end, tElapsed1 = toc(tStart1);
tStart2 = tic; for i=1:10000, ynext2 = local_state_space_iteration_k(yhat, epsilon, dr, M_, options_); end, tElapsed2 = toc(tStart2);