Fix `simult_.m` to call specialized simulation codes for order <= 3 instead of...

Fix `simult_.m` to call specialized simulation codes for order <= 3 instead of `k_order_simul` whenever possible

matlab/stochastic_solver/simult_.m

@@ -38,35 +38,28 @@ iter = size(ex_,1);
 endo_nbr = M_.endo_nbr;
 exo_nbr = M_.exo_nbr;
 
-y_ = zeros(size(y0,1),iter+M_.maximum_lag);
-y_(:,1) = y0;
-
 if options_.loglinear && ~options_.logged_steady_state
     k = get_all_variables_but_lagged_leaded_exogenous(M_);
     dr.ys(k)=log(dr.ys(k));
 end
 
-if ~options_.k_order_solver || (options_.k_order_solver && options_.pruning) %if k_order_pert is not used or if we do not use Dynare++ with k_order_pert
-    if iorder==1
-        y_(:,1) = y_(:,1)-dr.ys;
-    end
-end
-
-if options_.k_order_solver
+if (iorder > 3) || (iorder == 3 && ~options_.pruning)
     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'])
+        error(['The k_order_simul routine requires the specified approximation order to be '...
+               'consistent with the one used for computing the decision rules'])
     end
     y_ = k_order_simul(iorder,M_.nstatic,M_.npred,M_.nboth,M_.nfwrd,exo_nbr, ...
                        y0(dr.order_var,:),ex_',dr.ys(dr.order_var),dr, ...
                        options_.pruning);
     y_(dr.order_var,:) = y_;
 else
+    y_ = zeros(size(y0,1),iter+M_.maximum_lag);
+    y_(:,1) = y0;
     k2 = M_.nstatic+(1:M_.nspred);
     order_var = dr.order_var;
-
     switch iorder
       case 1
+        y_(:,1) = y_(:,1)-dr.ys;
         if isempty(dr.ghu)% For (linearized) deterministic models.
             for i = 2:iter+M_.maximum_lag
                 yhat = y_(order_var(k2),i-1);
@@ -163,7 +156,5 @@ else
             yhat2 = yhat2(ipred);
             yhat3 = yhat3(ipred);
         end
-      otherwise
-          error(['k_order_solver required for pruning at order = ' int2str(iorder)])
     end
 end

tests/stochastic_simulations/fs2000_k_order_simul.mod

@@ -65,13 +65,15 @@ end;
 
 steady;
 
-stoch_simul(order=2,nograph);
+stoch_simul(order=2,k_order_solver,nograph);
 
 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);
-Y2_mex = simult_(M_,options_,oo_.dr.ys,oo_.dr,ex,options_.order); 
+Y2_mex = k_order_simul(options_.order,M_.nstatic,M_.npred,M_.nboth,M_.nfwrd,M_.exo_nbr,   ...
+                       oo_.dr.ys(oo_.dr.order_var),ex',oo_.dr.ys(oo_.dr.order_var),oo_.dr,...
+                       options_.pruning);
+Y2_mex(oo_.dr.order_var,:) = Y2_mex;
 
 if max(abs(Y2_matlab(:)-Y2_mex(:)))>1e-8;
    error('2nd order: Matlab and mex simulation routines do not return similar results')
@@ -80,7 +82,10 @@ 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_temp = k_order_simul(options_.order,M_.nstatic,M_.npred,M_.nboth,M_.nfwrd, ...
+                            M_.exo_nbr,Y2_mexiter(oo_.dr.order_var,it),ex(it,:)', ...
+                            oo_.dr.ys(oo_.dr.order_var),oo_.dr,options_.pruning);
+    Y2_temp(oo_.dr.order_var,:) = Y2_temp;
     Y2_mexiter(:,it+1) = Y2_temp(:,2);    
 end
 
@@ -88,13 +93,18 @@ 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); 
-
+stoch_simul(order=3,k_order_solver,nograph,irf=0);
+Y3_mex = k_order_simul(options_.order,M_.nstatic,M_.npred,M_.nboth,M_.nfwrd,M_.exo_nbr,   ...
+                       oo_.dr.ys(oo_.dr.order_var),ex',oo_.dr.ys(oo_.dr.order_var),oo_.dr,...
+                       options_.pruning);
+Y3_mex(oo_.dr.order_var,:) = Y3_mex;
 Y3_mexiter = NaN(M_.endo_nbr,size(ex,1)+1);
 Y3_mexiter(:,1) = oo_.dr.ys;
 for it = 1:size(ex,1)
-    Y3_temp = simult_(M_,options_,Y3_mexiter(:,it),oo_.dr,ex(it,:),options_.order);
+    Y3_temp = k_order_simul(options_.order,M_.nstatic,M_.npred,M_.nboth,M_.nfwrd, ...
+                            M_.exo_nbr,Y3_mexiter(oo_.dr.order_var,it),ex(it,:)', ...
+                            oo_.dr.ys(oo_.dr.order_var),oo_.dr,options_.pruning);
+    Y3_temp(oo_.dr.order_var,:) = Y3_temp;
     Y3_mexiter(:,it+1) = Y3_temp(:,2);    
 end
 
@@ -102,23 +112,25 @@ 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);
+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); 
+Y2_mex = k_order_simul(options_.order,M_.nstatic,M_.npred,M_.nboth,M_.nfwrd,M_.exo_nbr,   ...
+                       oo_.dr.ys(oo_.dr.order_var),ex',oo_.dr.ys(oo_.dr.order_var),oo_.dr,...
+                       options_.pruning);
+Y2_mex(oo_.dr.order_var,:) = Y2_mex;
 
 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);
+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); 
+Y3_mex = k_order_simul(options_.order,M_.nstatic,M_.npred,M_.nboth,M_.nfwrd,M_.exo_nbr,   ...
+                       oo_.dr.ys(oo_.dr.order_var),ex',oo_.dr.ys(oo_.dr.order_var),oo_.dr,...
+                       options_.pruning);
+Y3_mex(oo_.dr.order_var,:) = Y3_mex;
 
 if max(abs(Y3_matlab(:)-Y3_mex(:)))>1e-8;
    error('3rd order with pruning: Matlab and mex simulation routines do not return similar results')