local_state_space_iteration_k MEX: fix bug in output

In its output, the MEX was returning values for all endogenous variables, but
it was used in a context where only the variables from oo_.dr.restrict_var_list
were expected (as is done with local_state_space_iteration_2 MEX).

This commit fixes this discrepancy, and also fixes the test that was checking
that both MEX are returning the same output.

Closes: #1768

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);
+
+        dr.eval(DecisionRule::emethod::horner, ynext_col_allvars, dyu);
+
+        ynext_col_allvars.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)};
-
-        dr.eval(DecisionRule::emethod::horner, ynext_col, dyu);
-
-        ynext_col.add(1.0, ys_reordered);
+        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();
 }

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

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.
+
+  This file must be run after first_spec.mod (both are based on the same model).
 */
 
-var y, k, a, h, b, c;
-varexo e, u;
+@#include "first_spec_common.inc"
 
-parameters beta, rho, alpha, delta, theta, psi, tau;
-
-alpha = 0.36;
-rho   = 0.95;
-tau   = 0.025;
-beta  = 0.99;
-delta = 0.025;
-psi   = 0;
-theta = 2.95;
-
-phi   = 0.1;
-
-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;
-end;
-
-initval;
-y = 1.08068253095672;
-c = 0.80359242014163;
-h = 0.29175631001732;
-k = 11.08360443260358;
-a = 0;
-b = 0;
-e = 0;
-u = 0;
-end;
+varobs q ca;
 
 shocks;
-var e; stderr 0.009;
-var u; stderr 0.009;
-var e, u = phi*0.009*0.009;
+var eeps = 0.04^2;
+var nnu = 0.03^2;
+var q = 0.01^2;
+var ca = 0.01^2;
 end;
 
-stoch_simul(order=2, irf=0, k_order_solver);
+// 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
+);
+
+stoch_simul(order=2, periods=200, 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);