diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst
index 2ebc3b4dc38f74dcf50853e070b3d96645a07c10..bebf95f34e481967c64533d75adbe0c3f6aa63b0 100644
--- a/doc/manual/source/the-model-file.rst
+++ b/doc/manual/source/the-model-file.rst
@@ -3011,27 +3011,27 @@ Finding the steady state with Dynare nonlinear solver
``12``
- Specialized version of ``2`` for models where all the equations
- have one endogenous variable on the left hand side and where
- each equation determines a different endogenous variable. Only
- expressions allowed on the left hand side are the natural
- logarithm of an endogenous variable, the first difference of an
- endogenous variable (with the ``diff`` operator), or the first
- difference of the logarithm of an endogenous variable.
- Univariate blocks are solved by evaluating the expression on the
- right hand side.
+ Computes a block decomposition and then applies a Newton-type
+ solver on those smaller blocks rather than on the full
+ nonlinear system. This is similar to ``2``, but is typically
+ more efficient. The block decomposition is done at the
+ preprocessor level, which brings two benefits: it identifies
+ blocks that can be evaluated rather than solved; and evulations
+ of the residual and Jacobian of the model are more efficient
+ because only the relevant elements are recomputed at every
+ iteration.
``14``
- Specialized version of ``4`` for models where all the equations
- have one endogenous variable on the left hand side and where
- each equation determines a different endogenous variable. Only
- expressions allowed on the left hand side are the natural
- logarithm of an endogenous variable, the first difference of an
- endogenous variable (with the ``diff`` operator), or the first
- difference of the logarithm of an endogenous variable..
- Univariate blocks are solved by evaluating the expression on the
- right hand side.
+ Computes a block decomposition and then applies a trust region
+ solver with autoscaling on those smaller blocks rather than on
+ the full nonlinear system. This is similar to ``4``, but is
+ typically more efficient. The block decomposition is done at
+ the preprocessor level, which brings two benefits: it
+ identifies blocks that can be evaluated rather than solved; and
+ evulations of the residual and Jacobian of the model are more
+ efficient because only the relevant elements are recomputed at
+ every iteration.
|br| Default value is ``4``.
diff --git a/matlab/backward/simul_backward_model.m b/matlab/backward/simul_backward_model.m
index 10333462401df897895e83213af902fdb520d11b..eb4f13088842dc7baca20d4439a4a54359e77ea2 100644
--- a/matlab/backward/simul_backward_model.m
+++ b/matlab/backward/simul_backward_model.m
@@ -47,10 +47,6 @@ if ~M_.maximum_lag
return
end
-if ismember(options_.solve_algo, [12,14]) && ~M_.possible_to_use_solve_algo_12_14
- error(M_.message_solve_algo_12_14)
-end
-
if nargin<3
Innovations = [];
else
@@ -100,4 +96,4 @@ if options_.linear
simulation = simul_backward_linear_model(initialconditions, samplesize, options_, M_, oo_, Innovations);
else
simulation = simul_backward_nonlinear_model(initialconditions, samplesize, options_, M_, oo_, Innovations);
-end
\ No newline at end of file
+end
diff --git a/matlab/backward/simul_backward_nonlinear_model_.m b/matlab/backward/simul_backward_nonlinear_model_.m
index 084ff04522de383bf975d8734f6138f4be39571a..c369d5d16f6d6411177a9c3b0acea9dbfad468a8 100644
--- a/matlab/backward/simul_backward_nonlinear_model_.m
+++ b/matlab/backward/simul_backward_nonlinear_model_.m
@@ -40,36 +40,64 @@ function [ysim, xsim] = simul_backward_nonlinear_model_(initialconditions, sampl
debug = false;
-model_dynamic_s = str2func('dynamic_backward_model_for_simulation');
-
if ~isempty(innovations)
DynareOutput.exo_simul(initialconditions.nobs+(1:samplesize),:) = innovations;
end
+if ismember(DynareOptions.solve_algo, [12,14])
+ [funcs, feedback_vars_idxs] = setup_time_recursive_block_simul(DynareModel);
+end
+
+function [r, J] = block_wrapper(z, feedback_vars_idx, func, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T)
+ % NB: do as few computations as possible inside this function, since it is
+ % called a very large number of times
+ y_dynamic(feedback_vars_idx) = z;
+ [~, ~, r, J] = feval(func, y_dynamic, x, DynareModel.params, DynareOutput.steady_state, ...
+ sparse_rowval, sparse_colval, sparse_colptr, T);
+end
+
% Simulations (call a Newton-like algorithm for each period).
for it = initialconditions.nobs+(1:samplesize)
if debug
dprintf('Period t = %s.', num2str(it-initialconditions.nobs));
end
y_ = DynareOutput.endo_simul(:,it-1);
- ylag = y_(iy1); % Set lagged variables.
y = y_; % A good guess for the initial conditions is the previous values for the endogenous variables.
try
if ismember(DynareOptions.solve_algo, [12,14])
- [DynareOutput.endo_simul(:,it), errorflag, ~, ~, errorcode] = dynare_solve(model_dynamic_s, y, ...
- DynareOptions.simul.maxit, DynareOptions.dynatol.f, DynareOptions.dynatol.x, ...
- DynareOptions, ...
- DynareModel.isloggedlhs, DynareModel.isauxdiffloggedrhs, DynareModel.endo_names, DynareModel.lhs, ...
- model_dynamic, ylag, DynareOutput.exo_simul, DynareModel.params, DynareOutput.steady_state, it);
+ x = DynareOutput.exo_simul(it,:);
+ T = NaN(DynareModel.block_structure.dyn_tmp_nbr);
+ y_dynamic = [y_; y; NaN(DynareModel.endo_nbr, 1)];
+ for blk = 1:length(DynareModel.block_structure.block)
+ sparse_rowval = DynareModel.block_structure.block(blk).g1_sparse_rowval;
+ sparse_colval = DynareModel.block_structure.block(blk).g1_sparse_colval;
+ sparse_colptr = DynareModel.block_structure.block(blk).g1_sparse_colptr;
+ if DynareModel.block_structure.block(blk).Simulation_Type ~= 1 % Not an evaluate forward block
+ [z, errorflag, ~, ~, errorcode] = dynare_solve(@block_wrapper, y_dynamic(feedback_vars_idxs{blk}), ...
+ DynareOptions.simul.maxit, DynareOptions.dynatol.f, ...
+ DynareOptions.dynatol.x, DynareOptions, ...
+ feedback_vars_idxs{blk}, funcs{blk}, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T);
+ if errorflag
+ error('Nonlinear solver routine failed with errorcode=%i in block %i and period %i.', errorcode, blk, it)
+ end
+ y_dynamic(feedback_vars_idxs{blk}) = z;
+ end
+ %% Compute endogenous if the block is of type evaluate or if there are recursive variables in a solve block.
+ %% Also update the temporary terms vector.
+ [y_dynamic, T] = feval(funcs{blk}, y_dynamic, x, DynareModel.params, ...
+ DynareOutput.steady_state, sparse_rowval, sparse_colval, ...
+ sparse_colptr, T);
+ end
+ DynareOutput.endo_simul(:,it) = y_dynamic(DynareModel.endo_nbr+(1:DynareModel.endo_nbr));
else
[DynareOutput.endo_simul(:,it), errorflag, ~, ~, errorcode] = ...
- dynare_solve(model_dynamic_s, y, ...
+ dynare_solve(@dynamic_backward_model_for_simulation, y, ...
DynareOptions.simul.maxit, DynareOptions.dynatol.f, DynareOptions.dynatol.x, ...
DynareOptions, ...
- model_dynamic, ylag, DynareOutput.exo_simul, DynareModel.params, DynareOutput.steady_state, it);
- end
- if errorflag
- error('Nonlinear solver routine failed with errorcode=%i in period %i.', errorcode, it)
+ model_dynamic, y_(iy1), DynareOutput.exo_simul, DynareModel.params, DynareOutput.steady_state, it);
+ if errorflag
+ error('Nonlinear solver routine failed with errorcode=%i in period %i.', errorcode, it)
+ end
end
catch
DynareOutput.endo_simul = DynareOutput.endo_simul(:, 1:it-1);
@@ -111,7 +139,7 @@ for it = initialconditions.nobs+(1:samplesize)
%
% Evaluate and check the residuals
%
- [r, J] = feval(model_dynamic_s, ytm, model_dynamic, ytm(iy1), DynareOutput.exo_simul, DynareModel.params, DynareOutput.steady_state, it);
+ [r, J] = feval(@dynamic_backward_model_for_simulation, ytm, model_dynamic, ytm(iy1), DynareOutput.exo_simul, DynareModel.params, DynareOutput.steady_state, it);
residuals_evaluating_to_nan = isnan(r);
residuals_evaluating_to_inf = isinf(r);
residuals_evaluating_to_complex = ~isreal(r);
@@ -142,6 +170,8 @@ end
ysim = DynareOutput.endo_simul(1:DynareModel.orig_endo_nbr,:);
xsim = DynareOutput.exo_simul;
+end
+
function display_names_of_problematic_equations(DynareModel, eqtags, TruthTable)
for i=1:DynareModel.orig_endo_nbr
if TruthTable(i)
@@ -152,4 +182,5 @@ for i=DynareModel.orig_endo_nbr+1:DynareModel.endo_nbr
if TruthTable(i)
dprintf(' - Auxiliary equation for %s', DynareModel.endo_names{i})
end
-end
\ No newline at end of file
+end
+end
diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index 17934f30ee560ff7966432599c6fec6d943b40bd..04266e61832b0e691fc8f43db081a52829133c92 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -58,21 +58,10 @@ else
in0 = nn;
end
-% Get first element of varargin if solve_algo ∈ {12,14} and rename varargin.
-if ismember(options.solve_algo, [12, 14])
- isloggedlhs = varargin{1};
- isauxdiffloggedrhs = varargin{2};
- endo_names = varargin{3};
- lhs = varargin{4};
- args = varargin(5:end);
-else
- args = varargin;
-end
-
% checking initial values
% TODO We should have an option to deactivate the randomization.
if jacobian_flag
- [fvec, fjac] = feval(f, x, args{:});
+ [fvec, fjac] = feval(f, x, varargin{:});
wrong_initial_guess_flag = false;
if ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:)))) || any(~isreal(fvec)) || any(~isreal(fjac(:)))
if ~ismember(options.solve_algo,[10,11]) && max(abs(fvec))< tolf
@@ -88,7 +77,7 @@ if jacobian_flag
while wrong_initial_guess_flag && tentative_number<=in0*10
tentative_number = tentative_number+1;
x(idx) = rand(in0, 1)*10;
- [fvec, fjac] = feval(f, x, args{:});
+ [fvec, fjac] = feval(f, x, varargin{:});
wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:))));
end
% If all previous attempts failed, try with real numbers.
@@ -96,7 +85,7 @@ if jacobian_flag
while wrong_initial_guess_flag && tentative_number<=in0*10
tentative_number = tentative_number+1;
x(idx) = randn(in0, 1)*10;
- [fvec, fjac] = feval(f, x, args{:});
+ [fvec, fjac] = feval(f, x, varargin{:});
wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:))));
end
% Last tentative, ff all previous attempts failed, try with negative numbers.
@@ -104,12 +93,12 @@ if jacobian_flag
while wrong_initial_guess_flag && tentative_number<=in0*10
tentative_number = tentative_number+1;
x(idx) = -rand(in0, 1)*10;
- [fvec, fjac] = feval(f, x, args{:});
+ [fvec, fjac] = feval(f, x, varargin{:});
wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:))));
end
end
else
- fvec = feval(f, x, args{:});
+ fvec = feval(f, x, varargin{:});
fjac = zeros(nn, nn);
if ~ismember(options.solve_algo,[10,11]) && max(abs(fvec)) < tolf
% return if initial value solves the problem except if a mixed complementarity problem is to be solved (complementarity conditions may not be satisfied)
@@ -125,7 +114,7 @@ else
while wrong_initial_guess_flag && tentative_number<=in0*10
tentative_number = tentative_number+1;
x(idx) = rand(in0, 1)*10;
- fvec = feval(f, x, args{:});
+ fvec = feval(f, x, varargin{:});
wrong_initial_guess_flag = ~all(isfinite(fvec));
end
% If all previous attempts failed, try with real numbers.
@@ -133,7 +122,7 @@ else
while wrong_initial_guess_flag && tentative_number<=in0*10
tentative_number = tentative_number+1;
x(idx) = randn(in0, 1)*10;
- fvec = feval(f, x, args{:});
+ fvec = feval(f, x, varargin{:});
wrong_initial_guess_flag = ~all(isfinite(fvec));
end
% Last tentative, ff all previous attempts failed, try with negative numbers.
@@ -141,7 +130,7 @@ else
while wrong_initial_guess_flag && tentative_number<=in0*10
tentative_number = tentative_number+1;
x(idx) = -rand(in0, 1)*10;
- fvec = feval(f, x, args{:});
+ fvec = feval(f, x, varargin{:});
wrong_initial_guess_flag = ~all(isfinite(fvec));
end
end
@@ -199,7 +188,7 @@ if options.solve_algo == 0
end
if ~isoctave
- [x, fvec, errorcode, ~, fjac] = fsolve(f, x, options4fsolve, args{:});
+ [x, fvec, errorcode, ~, fjac] = fsolve(f, x, options4fsolve, varargin{:});
else
% Under Octave, use a wrapper, since fsolve() does not have a 4th arg
if ischar(f)
@@ -207,7 +196,7 @@ if options.solve_algo == 0
else
f2 = f;
end
- [x, fvec, errorcode, ~, fjac] = fsolve(@(x) f2(x, args{:}), x, options4fsolve);
+ [x, fvec, errorcode, ~, fjac] = fsolve(@(x) f2(x, varargin{:}), x, options4fsolve);
end
if errorcode==1
errorflag = false;
@@ -220,91 +209,42 @@ if options.solve_algo == 0
else
errorflag = true;
end
-elseif options.solve_algo==1
- [x, errorflag, errorcode] = solve1(f, x, 1:nn, 1:nn, jacobian_flag, options.gstep, tolf, tolx, maxit, [], options.debug, args{:});
- [fvec, fjac] = feval(f, x, args{:});
+elseif ismember(options.solve_algo, [1, 12])
+ %% NB: It is the responsibility of the caller to deal with the block decomposition if solve_algo=12
+ [x, errorflag, errorcode] = solve1(f, x, 1:nn, 1:nn, jacobian_flag, options.gstep, tolf, tolx, maxit, [], options.debug, varargin{:});
+ [fvec, fjac] = feval(f, x, varargin{:});
elseif options.solve_algo==9
- [x, errorflag, errorcode] = trust_region(f, x, 1:nn, 1:nn, jacobian_flag, options.gstep, tolf, tolx, maxit, options.trust_region_initial_step_bound_factor, options.debug, args{:});
- [fvec, fjac] = feval(f, x, args{:});
-elseif ismember(options.solve_algo, [2, 12, 4])
- if ismember(options.solve_algo, [2, 12])
+ [x, errorflag, errorcode] = trust_region(f, x, 1:nn, 1:nn, jacobian_flag, options.gstep, tolf, tolx, maxit, options.trust_region_initial_step_bound_factor, options.debug, varargin{:});
+ [fvec, fjac] = feval(f, x, varargin{:});
+elseif ismember(options.solve_algo, [2, 4])
+ if options.solve_algo == 2
solver = @solve1;
else
solver = @trust_region;
end
- specializedunivariateblocks = options.solve_algo == 12;
if ~jacobian_flag
fjac = zeros(nn,nn) ;
dh = max(abs(x), options.gstep(1)*ones(nn,1))*eps^(1/3);
for j = 1:nn
xdh = x ;
xdh(j) = xdh(j)+dh(j) ;
- fjac(:,j) = (feval(f, xdh, args{:})-fvec)./dh(j) ;
+ fjac(:,j) = (feval(f, xdh, varargin{:})-fvec)./dh(j) ;
end
end
[j1,j2,r,s] = dmperm(fjac);
- JAC = abs(fjac(j1,j2))>0;
if options.debug
- disp(['DYNARE_SOLVE (solve_algo=2|4|12): number of blocks = ' num2str(length(r)-1)]);
+ disp(['DYNARE_SOLVE (solve_algo=2|4): number of blocks = ' num2str(length(r)-1)]);
end
- l = 0;
- fre = false;
for i=length(r)-1:-1:1
blocklength = r(i+1)-r(i);
if options.debug
- dprintf('DYNARE_SOLVE (solve_algo=2|4|12): solving block %u of size %u.', i, blocklength);
+ dprintf('DYNARE_SOLVE (solve_algo=2|4): solving block %u of size %u.', i, blocklength);
end
j = r(i):r(i+1)-1;
- if specializedunivariateblocks
- if options.debug
- dprintf('DYNARE_SOLVE (solve_algo=2|4|12): solving block %u by evaluating RHS.', i);
- end
- if isequal(blocklength, 1)
- if i<length(r)-1
- if fre || any(JAC(r(i), s(i)+(1:l)))
- % Reevaluation of the residuals is required because the current RHS depends on
- % variables that potentially have been updated previously.
- z = feval(f, x, args{:});
- l = 0;
- fre = false;
- end
- else
- % First iteration requires the evaluation of the residuals.
- z = feval(f, x, args{:});
- end
- l = l+1;
- if isequal(lhs{j1(j)}, endo_names{j2(j)}) || isequal(lhs{j1(j)}, sprintf('log(%s)', endo_names{j2(j)}))
- if isloggedlhs(j1(j))
- x(j2(j)) = exp(log(x(j2(j)))-z(j1(j)));
- else
- x(j2(j)) = x(j2(j))-z(j1(j));
- end
- else
- if options.debug
- dprintf('LHS variable is not determined by RHS expression (%u).', j1(j))
- dprintf('%s -> %s', lhs{j1(j)}, endo_names{j2(j)})
- end
- if ~isempty(regexp(lhs{j1(j)}, '\<AUX_DIFF_(\d*)\>', 'once'))
- if isauxdiffloggedrhs(j1(j))
- x(j2(j)) = exp(log(x(j2(j)))+z(j1(j)));
- else
- x(j2(j)) = x(j2(j))+z(j1(j));
- end
- else
- error('Algorithm solve_algo=%u cannot be used with this nonlinear problem.', options.solve_algo)
- end
- end
- continue
- end
- else
- if options.debug
- dprintf('DYNARE_SOLVE (solve_algo=2|4|12): solving block %u with trust_region routine.', i);
- end
- end
blockcolumns=s(i+1)-s(i);
if blockcolumns ~= blocklength
%non-square-block in DM; check whether initial value is solution
- [fval_check, fjac] = feval(f, x, args{:});
+ [fval_check, fjac] = feval(f, x, varargin{:});
if norm(fval_check(j1(j))) < tolf
errorflag = false;
errorcode = 0;
@@ -317,10 +257,9 @@ elseif ismember(options.solve_algo, [2, 12, 4])
options.gstep, ...
tolf, options.solve_tolx, maxit, ...
options.trust_region_initial_step_bound_factor, ...
- options.debug, args{:});
- fre = true;
+ options.debug, varargin{:});
else
- fprintf('\nDYNARE_SOLVE (solve_algo=2|4|12): the Dulmage-Mendelsohn decomposition returned a non-square block. This means that the Jacobian is singular. You may want to try another value for solve_algo.\n')
+ fprintf('\nDYNARE_SOLVE (solve_algo=2|4): the Dulmage-Mendelsohn decomposition returned a non-square block. This means that the Jacobian is singular. You may want to try another value for solve_algo.\n')
%overdetermined block
errorflag = true;
errorcode = 0;
@@ -329,31 +268,31 @@ elseif ismember(options.solve_algo, [2, 12, 4])
return
end
end
- fvec = feval(f, x, args{:});
+ fvec = feval(f, x, varargin{:});
if max(abs(fvec))>tolf
disp_verbose('Call solver on the full nonlinear problem.',options.verbosity)
[x, errorflag, errorcode] = solver(f, x, 1:nn, 1:nn, jacobian_flag, ...
options.gstep, tolf, options.solve_tolx, maxit, ...
options.trust_region_initial_step_bound_factor, ...
- options.debug, args{:});
+ options.debug, varargin{:});
end
- [fvec, fjac] = feval(f, x, args{:});
+ [fvec, fjac] = feval(f, x, varargin{:});
elseif options.solve_algo==3
if jacobian_flag
- [x, errorcode] = csolve(f, x, f, tolf, maxit, args{:});
+ [x, errorcode] = csolve(f, x, f, tolf, maxit, varargin{:});
else
- [x, errorcode] = csolve(f, x, [], tolf, maxit, args{:});
+ [x, errorcode] = csolve(f, x, [], tolf, maxit, varargin{:});
end
if errorcode==0
errorflag = false;
else
errorflag = true;
end
- [fvec, fjac] = feval(f, x, args{:});
+ [fvec, fjac] = feval(f, x, varargin{:});
elseif options.solve_algo==10
% LMMCP
olmmcp = options.lmmcp;
- [x, fvec, errorcode, ~, fjac] = lmmcp(f, x, olmmcp.lb, olmmcp.ub, olmmcp, args{:});
+ [x, fvec, errorcode, ~, fjac] = lmmcp(f, x, olmmcp.lb, olmmcp.ub, olmmcp, varargin{:});
eq_to_check=find(isfinite(olmmcp.lb) | isfinite(olmmcp.ub));
eq_to_ignore=eq_to_check(x(eq_to_check,:)<=olmmcp.lb(eq_to_check)+eps | x(eq_to_check,:)>=olmmcp.ub(eq_to_check)-eps);
fvec(eq_to_ignore)=0;
@@ -373,7 +312,7 @@ elseif options.solve_algo == 11
omcppath = options.mcppath;
global mcp_data
mcp_data.func = f;
- mcp_data.args = args;
+ mcp_data.args = varargin;
try
[x, fval, jac, mu] = pathmcp(x,omcppath.lb,omcppath.ub,'mcp_func',omcppath.A,omcppath.b,omcppath.t,omcppath.mu0);
catch
@@ -384,18 +323,15 @@ elseif options.solve_algo == 11
eq_to_ignore=eq_to_check(x(eq_to_check,:)<=omcppath.lb(eq_to_check)+eps | x(eq_to_check,:)>=omcppath.ub(eq_to_check)-eps);
fvec(eq_to_ignore)=0;
elseif ismember(options.solve_algo, [13, 14])
+ %% NB: It is the responsibility of the caller to deal with the block decomposition if solve_algo=14
if ~jacobian_flag
- error('DYNARE_SOLVE: option solve_algo=13|14 needs computed Jacobian')
- end
- auxstruct = struct();
- if options.solve_algo == 14
- auxstruct.lhs = lhs;
- auxstruct.endo_names = endo_names;
- auxstruct.isloggedlhs = isloggedlhs;
- auxstruct.isauxdiffloggedrhs = isauxdiffloggedrhs;
+ error('DYNARE_SOLVE: option solve_algo=13 needs computed Jacobian')
end
- [x, errorflag, errorcode] = block_trust_region(f, x, tolf, options.solve_tolx, maxit, options.trust_region_initial_step_bound_factor, options.debug, auxstruct, args{:});
- [fvec, fjac] = feval(f, x, args{:});
+ [x, errorflag, errorcode] = block_trust_region(f, x, tolf, options.solve_tolx, maxit, ...
+ options.trust_region_initial_step_bound_factor, ...
+ options.solve_algo == 13, ... % Only block-decompose with Dulmage-Mendelsohn for 13, not for 14
+ options.debug, varargin{:});
+ [fvec, fjac] = feval(f, x, varargin{:});
else
error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,9,10,11,12,13,14]')
end
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
index 65120b5f8af58a27e0631d690dfb07bcd9b4858b..7127000bf3e84b0ce5e024508c21b98f2297ec7e 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
@@ -50,6 +50,9 @@ if options_.linear && (isequal(options_.stack_solve_algo, 0) || isequal(options_
end
if options_.block
+ if M_.block_structure.time_recursive
+ error('Internal error: can''t perform stacked perfect foresight simulation with time-recursive block decomposition')
+ end
if options_.bytecode
try
oo_.endo_simul = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state,1, periods+2), periods);
diff --git a/matlab/perfect-foresight-models/setup_time_recursive_block_simul.m b/matlab/perfect-foresight-models/setup_time_recursive_block_simul.m
new file mode 100644
index 0000000000000000000000000000000000000000..d9b671b414c631c7c1798ba58e56799ec825d639
--- /dev/null
+++ b/matlab/perfect-foresight-models/setup_time_recursive_block_simul.m
@@ -0,0 +1,40 @@
+function [funcs, feedback_vars_idxs] = setup_time_recursive_block_simul(M_)
+%function [funcs, feedback_vars_idxs] = setup_time_recursive_block_simul(M_)
+%
+% For solve_algo={12,14}, precompute the function handles for per-block dynamic files
+% (it is surprisingly costly to recompute them within simulation loops).
+% Also precompute indices of feedback variables (also brings some performance gains).
+% By the way, do other sanity checks on block decomposition.
+
+% Copyright © 2022 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/>.
+
+if ~isfield(M_, 'block_structure')
+ error('Can''t use solve_algo=12 nor solve_algo=14, because the block decomposition of the dynamic model failed')
+end
+if ~M_.block_structure.time_recursive
+ error('Can''t use solve_algo=12 nor solve_algo=14, because the model is not purely backward/static/forward or you gave the ''block'' option to the ''model'' block')
+end
+
+funcs = cell(length(M_.block_structure.block), 1);
+feedback_vars_idxs = cell(length(M_.block_structure.block), 1);
+for blk = 1:length(M_.block_structure.block)
+ funcs{blk} = str2func(sprintf('%s.sparse.block.dynamic_%d', M_.fname, blk));
+ feedback_vars_idxs{blk} = M_.endo_nbr+M_.block_structure.block(blk).variable((M_.block_structure.block(blk).endo_nbr-M_.block_structure.block(blk).mfs+1):end); % Indices of feedback variables in the dynamic y vector (of size 3n)
+end
+
+end
diff --git a/matlab/perfect-foresight-models/sim1_purely_backward.m b/matlab/perfect-foresight-models/sim1_purely_backward.m
index ca0dcd011cb04debca708007a8ca1d27e197d130..c22a2d7ac198143a7bb234a3bbec15c8d2dc6463 100644
--- a/matlab/perfect-foresight-models/sim1_purely_backward.m
+++ b/matlab/perfect-foresight-models/sim1_purely_backward.m
@@ -20,38 +20,71 @@ function [endogenousvariables, info] = sim1_purely_backward(endogenousvariables,
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
ny0 = nnz(M.lead_lag_incidence(2,:)); % Number of variables at current period
-iyb = M.lead_lag_incidence(1,:)>0; % Logical vector (for lagged variables)
if ny0 ~= M.endo_nbr
error('All endogenous variables must appear at the current period!')
end
-if ismember(options.solve_algo, [12,14]) && ~M.possible_to_use_solve_algo_12_14
- error(M.message_solve_algo_12_14)
+if ismember(options.solve_algo, [12,14])
+ [funcs, feedback_vars_idxs] = setup_time_recursive_block_simul(M);
+else
+ iyb = M.lead_lag_incidence(1,:)>0; % Logical vector (for lagged variables)
+ dynamicmodel = str2func(sprintf('%s.%s', M.fname, 'dynamic'));
+ dynamicmodel_s = str2func('dynamic_backward_model_for_simulation');
end
-dynamicmodel = str2func(sprintf('%s.%s', M.fname, 'dynamic'));
-dynamicmodel_s = str2func('dynamic_backward_model_for_simulation');
+function [r, J] = block_wrapper(z, feedback_vars_idx, func, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T)
+ % NB: do as few computations as possible inside this function, since it is
+ % called a very large number of times
+ y_dynamic(feedback_vars_idx) = z;
+ [~, ~, r, J] = feval(func, y_dynamic, x, M.params, steadystate, ...
+ sparse_rowval, sparse_colval, sparse_colptr, T);
+end
info.status = true;
for it = M.maximum_lag + (1:options.periods)
y = endogenousvariables(:,it-1); % Values at previous period, also used as guess value for current period
- ylag = y(iyb);
if ismember(options.solve_algo, [12,14])
- [tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, y, ...
- options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
- options, M.isloggedlhs, M.isauxdiffloggedrhs, M.endo_names, M.lhs, ...
- dynamicmodel, ylag, exogenousvariables, M.params, steadystate, it);
+ x = exogenousvariables(it,:);
+ T = NaN(M.block_structure.dyn_tmp_nbr);
+ y_dynamic = [y; y; NaN(M.endo_nbr, 1)];
+ for blk = 1:length(M.block_structure.block)
+ sparse_rowval = M.block_structure.block(blk).g1_sparse_rowval;
+ sparse_colval = M.block_structure.block(blk).g1_sparse_colval;
+ sparse_colptr = M.block_structure.block(blk).g1_sparse_colptr;
+ if M.block_structure.block(blk).Simulation_Type ~= 1 % Not an evaluate forward block
+ [z, check, ~, ~, errorcode] = dynare_solve(@block_wrapper, y_dynamic(feedback_vars_idxs{blk}), ...
+ options.simul.maxit, options.dynatol.f, ...
+ options.dynatol.x, options, ...
+ feedback_vars_idxs{blk}, funcs{blk}, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T);
+ if check
+ info.status = false;
+ if options.debug
+ dprintf('sim1_purely_backward: Nonlinear solver routine failed with errorcode=%i in block %i and period %i.', errorcode, blk, it)
+ end
+ end
+ y_dynamic(feedback_vars_idxs{blk}) = z;
+ end
+ %% Compute endogenous if the block is of type evaluate or if there are recursive variables in a solve block.
+ %% Also update the temporary terms vector.
+ [y_dynamic, T] = feval(funcs{blk}, y_dynamic, x, M.params, ...
+ steadystate, sparse_rowval, sparse_colval, ...
+ sparse_colptr, T);
+ end
+ endogenousvariables(:,it) = y_dynamic(M.endo_nbr+(1:M.endo_nbr));
else
+ ylag = y(iyb);
[tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, y, ...
options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
options, dynamicmodel, ylag, exogenousvariables, M.params, steadystate, it);
+ if check
+ info.status = false;
+ dprintf('sim1_purely_backward: Nonlinear solver routine failed with errorcode=%i in period %i', errorcode, it)
+ break
+ end
+ endogenousvariables(:,it) = tmp;
end
- if check
- info.status = false;
- dprintf('sim1_purely_backward: Nonlinear solver routine failed with errorcode=%i in period %i', errorcode, it)
- break
- end
- endogenousvariables(:,it) = tmp;
-end
\ No newline at end of file
+end
+
+end
diff --git a/matlab/perfect-foresight-models/sim1_purely_forward.m b/matlab/perfect-foresight-models/sim1_purely_forward.m
index 715c544e817393fdbddefbd64da929bb811f5a53..a71715e2a60340f41b0a3a7beef41f2bec938201 100644
--- a/matlab/perfect-foresight-models/sim1_purely_forward.m
+++ b/matlab/perfect-foresight-models/sim1_purely_forward.m
@@ -19,34 +19,71 @@ function [endogenousvariables, info] = sim1_purely_forward(endogenousvariables,
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
ny0 = nnz(M.lead_lag_incidence(1,:)); % Number of variables at current period
-iyf = find(M.lead_lag_incidence(2,:)>0); % Indices of variables at next period
if ny0 ~= M.endo_nbr
error('All endogenous variables must appear at the current period!')
end
-dynamicmodel = str2func(sprintf('%s.%s', M.fname, 'dynamic'));
-dynamicmodel_s = str2func('dynamic_forward_model_for_simulation');
+if ismember(options.solve_algo, [12,14])
+ [funcs, feedback_vars_idxs] = setup_time_recursive_block_simul(M);
+else
+ iyf = find(M.lead_lag_incidence(2,:)>0); % Indices of variables at next period
+ dynamicmodel = str2func(sprintf('%s.%s', M.fname, 'dynamic'));
+ dynamicmodel_s = str2func('dynamic_forward_model_for_simulation');
+end
+
+function [r, J] = block_wrapper(z, feedback_vars_idx, func, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T)
+ % NB: do as few computations as possible inside this function, since it is
+ % called a very large number of times
+ y_dynamic(feedback_vars_idx) = z;
+ [~, ~, r, J] = feval(func, y_dynamic, x, M.params, steadystate, ...
+ sparse_rowval, sparse_colval, sparse_colptr, T);
+end
info.status = true;
for it = options.periods:-1:1
yf = endogenousvariables(:,it+1); % Values at next period, also used as guess value for current period
if ismember(options.solve_algo, [12,14])
- [tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, yf, ...
- options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
- options, M.isloggedlhs, M.isauxdiffloggedrhs, M.endo_names, M.lhs, ...
- dynamicmodel, yf(iyf), exogenousvariables, M.params, steadystate, it);
+ x = exogenousvariables(it,:);
+ T = NaN(M.block_structure.dyn_tmp_nbr);
+ y_dynamic = [NaN(M.endo_nbr, 1); yf; yf];
+ for blk = 1:length(M.block_structure.block)
+ sparse_rowval = M.block_structure.block(blk).g1_sparse_rowval;
+ sparse_colval = M.block_structure.block(blk).g1_sparse_colval;
+ sparse_colptr = M.block_structure.block(blk).g1_sparse_colptr;
+ if M.block_structure.block(blk).Simulation_Type ~= 2 % Not an evaluate backward block
+ [z, check, ~, ~, errorcode] = dynare_solve(@block_wrapper, y_dynamic(feedback_vars_idxs{blk}), ...
+ options.simul.maxit, options.dynatol.f, ...
+ options.dynatol.x, options, ...
+ feedback_vars_idxs{blk}, funcs{blk}, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T);
+ if check
+ info.status = false;
+ if options.debug
+ dprintf('sim1_purely_forward: Nonlinear solver routine failed with errorcode=%i in block %i and period %i.', errorcode, blk, it)
+ end
+ end
+ y_dynamic(feedback_vars_idxs{blk}) = z;
+ end
+ %% Compute endogenous if the block is of type evaluate or if there are recursive variables in a solve block.
+ %% Also update the temporary terms vector.
+ [y_dynamic, T] = feval(funcs{blk}, y_dynamic, x, M.params, ...
+ steadystate, sparse_rowval, sparse_colval, ...
+ sparse_colptr, T);
+ end
+ endogenousvariables(:,it) = y_dynamic(M.endo_nbr+(1:M.endo_nbr));
else
[tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, yf, ...
options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
options, ...
dynamicmodel, yf(iyf), exogenousvariables, M.params, steadystate, it);
+ if check
+ info.status = false;
+ dprintf('sim1_purely_forward: Nonlinear solver routine failed with errorcode=%i in period %i.', errorcode, it)
+ break
+ end
+ endogenousvariables(:,it) = tmp(1:M.endo_nbr);
end
- if check
- info.status = false;
- dprintf('sim1_purely_forward: Nonlinear solver routine failed with errorcode=%i in period %i.', errorcode, it)
- break
- end
- endogenousvariables(:,it) = tmp(1:M.endo_nbr);
-end
\ No newline at end of file
+end
+
+end
diff --git a/matlab/perfect-foresight-models/sim1_purely_static.m b/matlab/perfect-foresight-models/sim1_purely_static.m
index 8563d730d3319fbb2044a8075618d3fc7b678689..9e0c6623b7f3c8a55b4e07ae0b5ebee1ad1aae7c 100644
--- a/matlab/perfect-foresight-models/sim1_purely_static.m
+++ b/matlab/perfect-foresight-models/sim1_purely_static.m
@@ -23,12 +23,20 @@ if nnz(M.lead_lag_incidence(1,:)) ~= M.endo_nbr
error('All endogenous variables must appear at the current period!')
end
-if ismember(options.solve_algo, [12,14]) && ~M.possible_to_use_solve_algo_12_14
- error(M.message_solve_algo_12_14)
+if ismember(options.solve_algo, [12,14])
+ [funcs, feedback_vars_idxs] = setup_time_recursive_block_simul(M);
+else
+ dynamicmodel = str2func(sprintf('%s.%s', M.fname, 'dynamic'));
+ dynamicmodel_s = str2func('dynamic_static_model_for_simulation');
end
-dynamicmodel = str2func(sprintf('%s.%s', M.fname, 'dynamic'));
-dynamicmodel_s = str2func('dynamic_static_model_for_simulation');
+function [r, J] = block_wrapper(z, feedback_vars_idx, func, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T)
+ % NB: do as few computations as possible inside this function, since it is
+ % called a very large number of times
+ y_dynamic(feedback_vars_idx) = z;
+ [~, ~, r, J] = feval(func, y_dynamic, x, M.params, steadystate, ...
+ sparse_rowval, sparse_colval, sparse_colptr, T);
+end
info.status = true;
@@ -36,21 +44,46 @@ y = endogenousvariables(:,1);
for it = 1:options.periods
if ismember(options.solve_algo, [12,14])
- [tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, y, ...
- options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
- options, M.isloggedlhs, M.isauxdiffloggedrhs, M.endo_names, M.lhs, ...
- dynamicmodel, exogenousvariables, M.params, steadystate, it);
+ x = exogenousvariables(it,:);
+ T = NaN(M.block_structure.dyn_tmp_nbr);
+ y_dynamic = [NaN(M.endo_nbr, 1); y; NaN(M.endo_nbr, 1)];
+ for blk = 1:length(M.block_structure.block)
+ sparse_rowval = M.block_structure.block(blk).g1_sparse_rowval;
+ sparse_colval = M.block_structure.block(blk).g1_sparse_colval;
+ sparse_colptr = M.block_structure.block(blk).g1_sparse_colptr;
+ if M.block_structure.block(blk).Simulation_Type ~= 1 % Not an evaluate forward block
+ [z, check, ~, ~, errorcode] = dynare_solve(@block_wrapper, y_dynamic(feedback_vars_idxs{blk}), ...
+ options.simul.maxit, options.dynatol.f, ...
+ options.dynatol.x, options, ...
+ feedback_vars_idxs{blk}, funcs{blk}, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T);
+ if check
+ info.status = false;
+ if options.debug
+ dprintf('sim1_purely_static: Nonlinear solver routine failed with errorcode=%i in block %i and period %i.', errorcode, blk, it)
+ end
+ end
+ y_dynamic(feedback_vars_idxs{blk}) = z;
+ end
+ %% Compute endogenous if the block is of type evaluate or if there are recursive variables in a solve block.
+ %% Also update the temporary terms vector.
+ [y_dynamic, T] = feval(funcs{blk}, y_dynamic, x, M.params, ...
+ steadystate, sparse_rowval, sparse_colval, ...
+ sparse_colptr, T);
+ end
+ endogenousvariables(:,it) = y_dynamic(M.endo_nbr+(1:M.endo_nbr));
else
[tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, y, ...
options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
options, dynamicmodel, exogenousvariables, M.params, steadystate, it);
- end
- if check
- info.status = false;
- if options.debug
- dprintf('sim1_purely_static: Nonlinear solver routine failed with errorcode=%i in period %i.', errorcode, it)
+ if check
+ info.status = false;
+ if options.debug
+ dprintf('sim1_purely_static: Nonlinear solver routine failed with errorcode=%i in period %i.', errorcode, it)
+ end
end
+ endogenousvariables(:,it) = tmp;
end
- endogenousvariables(:,it) = tmp;
y = endogenousvariables(:,it);
-end
\ No newline at end of file
+end
+
+end
diff --git a/matlab/setup_solvers.m b/matlab/setup_solvers.m
deleted file mode 100644
index 13aac7b5326fdadd1e1bba2c50298b0f0cc4578d..0000000000000000000000000000000000000000
--- a/matlab/setup_solvers.m
+++ /dev/null
@@ -1,97 +0,0 @@
-function Model = setup_solvers(Model)
-
-% Setup solve_algo={12,14} by identifying equations with a log on the left hand side.
-%
-% INPUTS
-% - Model [struct] Model description, aka M_.
-% - Options [struct] Dynare's options, aka options_.
-%
-% OUTPUTS
-% - Model [struct] Updated model description.
-
-% Copyright © 2020 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/>.
-
-cannot_use_solve_algo_12_14 = false;
-
-try
- json = loadjson_(sprintf('%s/model/json/modfile.json', Model.fname));
-catch
- cannot_use_solve_algo_12_14 = true;
- message = 'Algorithms solve_algo={12,14} require json output of the model (use json=compute option)';
-end
-
-if ~cannot_use_solve_algo_12_14
-
- lhs = cell(length(json.model),1);
- isauxdiffloggedrhs = false(length(json.model), 1);
-
- for i = 1:length(json.model)
- if length(json.model)>1
- lhs{i} = json.model{i}.lhs;
- else
- lhs{i} = json.model.lhs;
- end
- if isempty(regexp(lhs{i}, '^\w+$|^log\(\w+\)$'))
- cannot_use_solve_algo_12_14 = true;
- message = sprintf('With solve_algo={12,14}, each equation must have on the left hand side a single variable or logged variable (equation %d does not satisfy this condition).', i);
- break
- end
- if length(json.model)>1
- rhs = json.model{i}.rhs;
- else
- rhs = json.model.lhs;
- end
- if i>Model.orig_endo_nbr && ~isempty(regexp(lhs{i}, '\<AUX_DIFF_(\d*)\>', 'once')) && ismember(lhs{i}, lhs(1:i-1)) && ...
- ~isempty(regexp(rhs, 'log\(\w*\)-log\(\w*\(-1\)\)', 'once'))
- isauxdiffloggedrhs(i) = true;
- end
- end
-
-end
-
-if ~cannot_use_solve_algo_12_14
-
- islog = @(x) ~isempty(regexp(x, 'log\(\w*\)', 'once'));
-
- lhs0 = lhs;
- for i=1:length(json.model)
- if islog(lhs{i})
- lhs0{i} = strrep(strrep(lhs{i}, 'log(', ''), ')', '');
- end
- end
-
- if ~isequal(length(unique(lhs0(1:Model.orig_endo_nbr))), length(lhs0(1:Model.orig_endo_nbr)))
- cannot_use_solve_algo_12_14 = true;
- message = sprintf('With solve_algo={12,14}, each equation must determine a different endogenous variable.')
- end
-
-end
-
-if cannot_use_solve_algo_12_14
- Model.isloggedlhs = {};
- Model.lhs = {};
- Model.isauxdiffloggedrhs = [];
- Model.possible_to_use_solve_algo_12_14 = false;
- Model.message_solve_algo_12_14 = message;
-else
- Model.isloggedlhs = cellfun(islog, lhs);
- Model.lhs = lhs;
- Model.isauxdiffloggedrhs = isauxdiffloggedrhs;
- Model.possible_to_use_solve_algo_12_14 = true;
- Model.message_solve_algo_12_14 = '';
-end
\ No newline at end of file
diff --git a/mex/sources/block_trust_region/mexFunction.f08 b/mex/sources/block_trust_region/mexFunction.f08
index 59da65547f86acda02deb9f9f241271b761e6330..11f9fccf209d330a87fc8c30fe9e9b725e6d6b76 100644
--- a/mex/sources/block_trust_region/mexFunction.f08
+++ b/mex/sources/block_trust_region/mexFunction.f08
@@ -35,16 +35,11 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
integer :: maxiter
real(real64), dimension(:), allocatable :: fvec
real(real64), dimension(:,:), allocatable :: fjac
- logical :: debug, specializedunivariateblocks
+ logical :: debug, block_decompose
character(len=80) :: debug_msg
- logical(mxLogical), dimension(:), pointer :: isloggedlhs => null(), &
- isauxdiffloggedrhs => null()
- type(c_ptr) :: endo_names, lhs
- logical :: fre ! True if the last block has been solved (i.e. not evaluated), so that residuals must be updated
- integer, dimension(:), allocatable :: evaled_cols ! If fre=.false., lists the columns that have been evaluated so far without updating the residuals
- if (nrhs < 4 .or. nlhs /= 3) then
- call mexErrMsgTxt("Must have at least 7 inputs and exactly 3 outputs")
+ if (nrhs < 8 .or. nlhs /= 3) then
+ call mexErrMsgTxt("Must have at least 8 inputs and exactly 3 outputs")
end if
if (.not. ((mxIsChar(prhs(1)) .and. mxGetM(prhs(1)) == 1) .or. mxIsClass(prhs(1), "function_handle"))) then
@@ -73,178 +68,96 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
end if
if (.not. (mxIsLogicalScalar(prhs(7)))) then
- call mexErrMsgTxt("Seventh argument (debug) should be a logical scalar")
+ call mexErrMsgTxt("Seventh argument (block_decompose) should be a logical scalar")
end if
- if (.not. (mxIsStruct(prhs(8)) .and. &
- (mxGetNumberOfFields(prhs(8)) == 0 .or. mxGetNumberOfFields(prhs(8)) == 4))) then
- call mexErrMsgTxt("Eighth argument should be a struct with either 0 or 4 fields")
+ if (.not. (mxIsLogicalScalar(prhs(8)))) then
+ call mexErrMsgTxt("Eigth argument (debug) should be a logical scalar")
end if
- specializedunivariateblocks = (mxGetNumberOfFields(prhs(8)) == 4)
-
func => prhs(1)
tolf = mxGetScalar(prhs(3))
tolx = mxGetScalar(prhs(4))
maxiter = int(mxGetScalar(prhs(5)))
factor = mxGetScalar(prhs(6))
- debug = mxGetScalar(prhs(7)) == 1._c_double
- extra_args => prhs(9:nrhs) ! Extra arguments to func are in argument 9 and subsequent ones
+ block_decompose = mxGetScalar(prhs(7)) == 1._c_double
+ debug = mxGetScalar(prhs(8)) == 1._c_double
+ extra_args => prhs(9:nrhs) ! Extra arguments to func are in argument 8 and subsequent ones
associate (x_mat => mxGetPr(prhs(2)))
allocate(x(size(x_mat)))
x = x_mat
end associate
- if (specializedunivariateblocks) then
- block
- type(c_ptr) :: tmp
- tmp = mxGetField(prhs(8), 1_mwIndex, "isloggedlhs")
- if (.not. (c_associated(tmp) .and. mxIsLogical(tmp) .and. mxGetNumberOfElements(tmp) == size(x))) then
- call mexErrMsgTxt("Seventh argument must have a 'isloggedlhs' field of type logical, of same size as second argument")
- end if
- isloggedlhs => mxGetLogicals(tmp)
-
- tmp = mxGetField(prhs(8), 1_mwIndex, "isauxdiffloggedrhs")
- if (.not. (c_associated(tmp) .and. mxIsLogical(tmp) .and. mxGetNumberOfElements(tmp) == size(x))) then
- call mexErrMsgTxt("Seventh argument must have a 'isauxdiffloggedrhs' field of type &
- &logical, of same size as second argument")
- end if
- isauxdiffloggedrhs => mxGetLogicals(tmp)
-
- lhs = mxGetField(prhs(8), 1_mwIndex, "lhs")
- if (.not. (c_associated(lhs) .and. mxIsCell(lhs) .and. mxGetNumberOfElements(lhs) == size(x))) then
- call mexErrMsgTxt("Seventh argument must have a 'lhs' field of type cell, of same size as second argument")
- end if
-
- endo_names = mxGetField(prhs(8), 1_mwIndex, "endo_names")
- if (.not. (c_associated(endo_names) .and. mxIsCell(endo_names) .and. mxGetNumberOfElements(endo_names) == size(x))) then
- call mexErrMsgTxt("Seventh argument must have a 'endo_names' field of type cell, of same size as second argument")
- end if
- end block
-
- allocate(evaled_cols(0))
- fre = .false.
- end if
allocate(fvec(size(x)))
allocate(fjac(size(x), size(x)))
- ! Compute block decomposition
- nullify(x_indices, f_indices, x_all)
- call matlab_fcn(x, fvec, fjac)
- call dm_blocks(fjac, blocks)
-
- if (debug) then
- write (debug_msg, "('DYNARE_SOLVE (solve_algo=13|14): number of blocks = ', i0)") size(blocks)
- call mexPrintf_trim_newline(debug_msg)
- end if
+ if (block_decompose) then
+ ! Compute block decomposition
+ nullify(x_indices, f_indices, x_all)
+ call matlab_fcn(x, fvec, fjac)
+ call dm_blocks(fjac, blocks)
- ! Solve the system, starting from bottom-rightmost block
- do i = size(blocks),1,-1
if (debug) then
- write (debug_msg, "('DYNARE_SOLVE (solve_algo=13|14): solving block ', i0, ' of size ', i0)") &
- i, size(blocks(i)%col_indices)
+ write (debug_msg, "('DYNARE_SOLVE (solve_algo=13): number of blocks = ', i0)") size(blocks)
call mexPrintf_trim_newline(debug_msg)
end if
- if (specializedunivariateblocks .and. size(blocks(i)%col_indices) == 1) then
+ ! Solve the system, starting from bottom-rightmost block
+ do i = size(blocks),1,-1
if (debug) then
- write (debug_msg, "('DYNARE_SOLVE (solve_algo=13|14): solving block ', i0, ' by evaluating RHS')") i
+ write (debug_msg, "('DYNARE_SOLVE (solve_algo=13): solving block ', &
+ &i0, ' of size ', i0)") &
+ i, size(blocks(i)%col_indices)
call mexPrintf_trim_newline(debug_msg)
end if
- associate (eq => blocks(i)%row_indices(1), var => blocks(i)%col_indices(1))
- if (fre .or. any(abs(fjac(eq, evaled_cols)) > 0._real64)) then
- ! Reevaluation of the residuals is required because the current RHS depends on
- ! variables that potentially have been updated previously.
- nullify(x_indices, f_indices, x_all)
- call matlab_fcn(x, fvec)
- deallocate(evaled_cols) ! This shouldn’t be necessary, but it crashes otherwise with gfortran 8
- allocate(evaled_cols(0))
- fre = .false.
- end if
- evaled_cols = [ evaled_cols, var]
- block
- ! An associate() construct for lhs_eq and endo_name_var makes the
- ! code crash (with double free) using gfortran 8. Hence use a block
- character(kind=c_char, len=:), allocatable :: lhs_eq, endo_name_var
- lhs_eq = mxArrayToString(mxGetCell(lhs, int(eq, mwIndex)))
- endo_name_var = mxArrayToString(mxGetCell(endo_names, int(var, mwIndex)))
- if (lhs_eq == endo_name_var .or. lhs_eq == "log(" // endo_name_var // ")") then
- if (isloggedlhs(eq)) then
- x(var) = exp(log(x(var)) - fvec(eq))
- else
- x(var) = x(var) - fvec(eq)
- end if
- else
- if (debug) then
- write (debug_msg, "('LHS variable is not determined by RHS expression (', i0, ')')") eq
- call mexPrintf_trim_newline(debug_msg)
- write (debug_msg, "(a, ' -> ', a)") lhs_eq, endo_name_var
- call mexPrintf_trim_newline(debug_msg)
- end if
- if (lhs_eq(1:9) == "AUX_DIFF_" .or. lhs_eq(1:13) == "log(AUX_DIFF_") then
- if (isauxdiffloggedrhs(eq)) then
- x(var) = exp(log(x(var)) + fvec(eq))
- else
- x(var) = x(var) + fvec(eq)
- end if
- else
- call mexErrMsgTxt("Algorithm solve_algo=14 cannot be used with this nonlinear problem")
- end if
- end if
- end block
- end associate
- cycle
- else
- if (debug) then
- write (debug_msg, "('DYNARE_SOLVE (solve_algo=13|14): solving block ', i0, ' with trust region routine')") i
- call mexPrintf_trim_newline(debug_msg)
+
+ if (size(blocks(i)%col_indices) /= size(blocks(i)%row_indices)) then
+ ! Non-square block in DM decomposition
+ ! Before erroring out, check whether we are not already at the solution for this block
+ ! See also #1851
+ if (norm2(fvec(blocks(i)%row_indices)) < tolf) then
+ cycle
+ else
+ call mexErrMsgTxt("DYNARE_SOLVE (solve_algo=13): the Dulmage-Mendelsohn &
+ &decomposition returned a non-square block. This means that the &
+ &Jacobian is singular. You may want to try another value for solve_algo.")
+ end if
end if
- end if
- if (size(blocks(i)%col_indices) /= size(blocks(i)%row_indices)) then
- ! Non-square block in DM decomposition
- ! Before erroring out, check whether we are not already at the solution for this block
- ! See also #1851
- if (norm2(fvec(blocks(i)%row_indices)) < tolf) then
- cycle
- else
- call mexErrMsgTxt("DYNARE_SOLVE (solve_algo=13|14): the Dulmage-Mendelsohn &
- &decomposition returned a non-square block. This means that the &
- &Jacobian is singular. You may want to try another value for solve_algo.")
+ block
+ real(real64), dimension(size(blocks(i)%col_indices)) :: x_block
+ x_indices => blocks(i)%col_indices
+ f_indices => blocks(i)%row_indices
+ x_all => x
+ x_block = x(x_indices)
+ call trust_region_solve(x_block, matlab_fcn, info, tolx, tolf, maxiter, factor)
+ x(x_indices) = x_block
+ end block
+ end do
+
+ ! Verify that we have a solution
+ ! Note that here we use the ∞-norm, while trust region uses 2-norm; otherwise
+ ! this check would almost always fail (because the 2-norm of the full fvec is
+ ! larger than the 2-norm of its sub-vectors)
+ ! If the check fails, this normally means that the block decomposition was
+ ! incorrect (because some element of the Jacobian was numerically zero at the
+ ! guess value while not being symbolically zero)
+ nullify(x_indices, f_indices, x_all)
+ call matlab_fcn(x, fvec)
+ if (maxval(abs(fvec)) > tolf) then
+ if (debug) &
+ call mexPrintf_trim_newline("DYNARE_SOLVE (solve_algo=13): residuals&
+ & still too large, solving for the whole model")
+ call trust_region_solve(x, matlab_fcn, info, tolx, tolf, maxiter, factor)
+ else
+ if (size(blocks) > 1) then
+ ! Note that the value of info may be different across blocks
+ info = 1
end if
end if
-
- block
- real(real64), dimension(size(blocks(i)%col_indices)) :: x_block
- x_indices => blocks(i)%col_indices
- f_indices => blocks(i)%row_indices
- x_all => x
- x_block = x(x_indices)
- call trust_region_solve(x_block, matlab_fcn, info, tolx, tolf, maxiter, factor)
- x(x_indices) = x_block
- end block
-
- fre = .true.
- end do
-
- ! Verify that we have a solution
- ! Note that here we use the ∞-norm, while trust region uses 2-norm; otherwise
- ! this check would almost always fail (because the 2-norm of the full fvec is
- ! larger than the 2-norm of its sub-vectors)
- ! If the check fails, this normally means that the block decomposition was
- ! incorrect (because some element of the Jacobian was numerically zero at the
- ! guess value while not being symbolically zero)
- nullify(x_indices, f_indices, x_all)
- call matlab_fcn(x, fvec)
- if (maxval(abs(fvec)) > tolf) then
- if (debug) &
- call mexPrintf_trim_newline("DYNARE_SOLVE (solve_algo=13|14): residuals still too large, solving for the whole model")
+ else ! No block decomposition
+ nullify(x_indices, f_indices, x_all)
call trust_region_solve(x, matlab_fcn, info, tolx, tolf, maxiter, factor)
- else
- if (size(blocks) > 1) then
- ! Note that the value of info may be different across blocks
- info = 1
- end if
end if
plhs(1) = mxCreateDoubleMatrix(int(size(x, 1), mwSize), 1_mwSize, mxREAL)
diff --git a/preprocessor b/preprocessor
index f725c534ef1344813421c2404de2cb4f3774d113..c48248fc0d65ed70fb13dd508ed0254f113cbcb6 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit f725c534ef1344813421c2404de2cb4f3774d113
+Subproject commit c48248fc0d65ed70fb13dd508ed0254f113cbcb6
diff --git a/tests/Makefile.am b/tests/Makefile.am
index b8b04355370eba4f5b7644e400a3b89f22ba69ea..e0f8aead52e2a8b3573ee6a1da6bf0a48fc43d16 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -498,7 +498,19 @@ MODFILES = \
log_transform/example1.mod \
log_transform/ramst.mod \
log_transform/fs2000_nonstationary.mod \
- log_transform/nk_ramsey.mod
+ log_transform/nk_ramsey.mod \
+ solve_algo_12_14/simul_backward_reference.mod \
+ solve_algo_12_14/simul_backward_12.mod \
+ solve_algo_12_14/simul_backward_14.mod \
+ solve_algo_12_14/purely_backward_reference.mod \
+ solve_algo_12_14/purely_backward_12.mod \
+ solve_algo_12_14/purely_backward_14.mod \
+ solve_algo_12_14/purely_static_reference.mod \
+ solve_algo_12_14/purely_static_12.mod \
+ solve_algo_12_14/purely_static_14.mod \
+ solve_algo_12_14/purely_forward_reference.mod \
+ solve_algo_12_14/purely_forward_12.mod \
+ solve_algo_12_14/purely_forward_14.mod
ECB_MODFILES = \
var-expectations/1/example1.mod \
@@ -994,6 +1006,27 @@ model-inversion/nk-2/z_check_inversion.o.trs: model-inversion/nk-2/invert.o.trs
ep/rbcii_MCP.m.trs: ep/rbcii.m.trs
ep/rbcii_MCP.o.trs: ep/rbcii.o.trs
+solve_algo_12_14/simul_backward_12.m.trs: solve_algo_12_14/simul_backward_reference.m.trs
+solve_algo_12_14/simul_backward_12.o.trs: solve_algo_12_14/simul_backward_reference.o.trs
+solve_algo_12_14/simul_backward_14.m.trs: solve_algo_12_14/simul_backward_reference.m.trs
+solve_algo_12_14/simul_backward_14.o.trs: solve_algo_12_14/simul_backward_reference.o.trs
+
+solve_algo_12_14/purely_backward_12.m.trs: solve_algo_12_14/purely_backward_reference.m.trs
+solve_algo_12_14/purely_backward_12.o.trs: solve_algo_12_14/purely_backward_reference.o.trs
+solve_algo_12_14/purely_backward_14.m.trs: solve_algo_12_14/purely_backward_reference.m.trs
+solve_algo_12_14/purely_backward_14.o.trs: solve_algo_12_14/purely_backward_reference.o.trs
+
+solve_algo_12_14/purely_static_12.m.trs: solve_algo_12_14/purely_static_reference.m.trs
+solve_algo_12_14/purely_static_12.o.trs: solve_algo_12_14/purely_static_reference.o.trs
+solve_algo_12_14/purely_static_14.m.trs: solve_algo_12_14/purely_static_reference.m.trs
+solve_algo_12_14/purely_static_14.o.trs: solve_algo_12_14/purely_static_reference.o.trs
+
+solve_algo_12_14/purely_forward_12.m.trs: solve_algo_12_14/purely_forward_reference.m.trs
+solve_algo_12_14/purely_forward_12.o.trs: solve_algo_12_14/purely_forward_reference.o.trs
+solve_algo_12_14/purely_forward_14.m.trs: solve_algo_12_14/purely_forward_reference.m.trs
+solve_algo_12_14/purely_forward_14.o.trs: solve_algo_12_14/purely_forward_reference.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)))
@@ -1210,6 +1243,11 @@ model-inversion: m/model-inversion o/model-inversion
m/model-inversion: $(patsubst %.mod, %.m.trs, $(filter model-inversion/%.mod, $(MODFILES)))
o/model-inversion: $(patsubst %.mod, %.o.trs, $(filter model-inversion/%.mod, $(MODFILES)))
+solve_algo_12_14: m/solve_algo_12_14 o/solve_algo_12_14
+m/solve_algo_12_14: $(patsubst %.mod, %.m.trs, $(filter solve_algo_12_14/%.mod, $(MODFILES)))
+o/solve_algo_12_14: $(patsubst %.mod, %.o.trs, $(filter solve_algo_12_14/%.mod, $(MODFILES)))
+
+
# ECB files
M_ECB_TRS_FILES = $(patsubst %.mod, %.m.trs, $(ECB_MODFILES))
@@ -1448,7 +1486,12 @@ EXTRA_DIST = \
solver-test-functions/variablydimensioned.m \
solver-test-functions/watson.m \
solver-test-functions/wood.m \
- ramst_normcdf_and_friends.inc
+ ramst_normcdf_and_friends.inc \
+ solve_algo_12_14/backward_model.inc \
+ solve_algo_12_14/simul_backward_common.inc \
+ solve_algo_12_14/purely_backward_common.inc \
+ solve_algo_12_14/purely_static_common.inc \
+ solve_algo_12_14/purely_forward_common.inc
if ENABLE_MATLAB
check-local: check-matlab
@@ -1639,6 +1682,8 @@ clean-local:
rm -rf tests/pac/var-12/toto
+ rm -f solve_algo_12_14/simul_backward_ref.mat
+
find . -name "*.tex" -type f -delete
find . -name "*.aux" -type f -delete
find . -name "*.log" -type f -delete
diff --git a/tests/nonlinearsolvers.m b/tests/nonlinearsolvers.m
index aab3e2a7a4c5c7528ed9f50fb288114a535f8163..b8e07f3c4a089466e317b1816262fdd131cef1f1 100644
--- a/tests/nonlinearsolvers.m
+++ b/tests/nonlinearsolvers.m
@@ -31,8 +31,6 @@ tolx = 1e-6;
maxit = 50;
factor = 10;
-auxstruct = struct();
-
% List of function handles
objfun = { @rosenbrock,
@powell1,
@@ -73,7 +71,7 @@ for i=1:length(objfun)
x = objfun{i}();
end
try
- [x, errorflag, exitflag] = block_trust_region(objfun{i}, x, tolf, tolx, maxit, factor, false, auxstruct);
+ [x, errorflag, exitflag] = block_trust_region(objfun{i}, x, tolf, tolx, maxit, factor, true, false);
if isequal(func2str(objfun{i}), 'powell2')
if ~errorflag
testFailed = testFailed+1;
diff --git a/tests/solve_algo_12_14/backward_model.inc b/tests/solve_algo_12_14/backward_model.inc
new file mode 100644
index 0000000000000000000000000000000000000000..43c097c68f7b314ac17f48b5a315dc99d04a24a3
--- /dev/null
+++ b/tests/solve_algo_12_14/backward_model.inc
@@ -0,0 +1,28 @@
+/* A simple purely backward model, used for testing both “simul_backward” and
+ perfect foresight simulations, consisting of:
+ - a VAR(3) with one variable in level, one in log and one log-differentiated;
+ - two other variables that depend on the contemporaneous values of the VAR(3),
+ and which together for a simultaneous block.
+ Those five equations are repeated an arbitrary number of times (n). */
+
+@#define n = 10
+
+@#for i in 1:n
+var x@{i}, y@{i}, z@{i}, t@{i}, s@{i};
+varexo e_x@{i}, e_y@{i}, e_z@{i}, e_t@{i}, e_s@{i};
+@#endfor
+
+model(use_dll);
+@#for i in 1:n
+x@{i} = 0.1*x@{i}(-1) + 0.2*log(y@{i}(-1)) + 0.3*diff(log(z@{i}(-1))) + e_x@{i};
+[name = 'y']
+log(y@{i}) = 0.4*x@{i}(-1) + 0.5*log(y@{i}(-1)) + 0.6*diff(log(z@{i}(-1))) + e_y@{i};
+[name = 'z']
+/* NB: we choose 0.5 as steady state for z, since initial value is drawn from
+ [0,1] (in the “simul_backward” case) */
+diff(log(z@{i})) = -0.8*(log(z@{i}(-1)) - log(0.5)) + 0.1*x@{i}(-1) + 0.2*log(y@{i}(-1)) - 0.3*diff(log(z@{i}(-1))) + e_z@{i};
+
+t@{i} = x@{i} + 2*log(y@{i}) + 3*diff(log(z@{i})) + s@{i} + e_t@{i};
+s@{i} = 4*x@{i}(-1) + 3*t@{i} + e_s@{i};
+@#endfor
+end;
diff --git a/tests/solve_algo_12_14/purely_backward_12.mod b/tests/solve_algo_12_14/purely_backward_12.mod
new file mode 100644
index 0000000000000000000000000000000000000000..ffcfbc9ed75a8704327381ca43cde488af965a11
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_backward_12.mod
@@ -0,0 +1,13 @@
+// Check the correctedness of perfect foresight simulation of a purely backward model with “solve_algo=12”
+
+@#include "purely_backward_common.inc"
+
+perfect_foresight_solver(solve_algo = 12);
+
+ref = load('purely_backward_reference/Output/purely_backward_reference_results.mat');
+
+if max(max(abs(oo_.endo_simul - ref.oo_.endo_simul))) > 5e-5
+ error('Incorrect results for perfect foresight with solve_algo=12')
+end
+
+
diff --git a/tests/solve_algo_12_14/purely_backward_14.mod b/tests/solve_algo_12_14/purely_backward_14.mod
new file mode 100644
index 0000000000000000000000000000000000000000..238c188552d3ad85c772766362d2b90f9e3ff90c
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_backward_14.mod
@@ -0,0 +1,13 @@
+// Check the correctedness of perfect foresight simulation of a purely backward model with “solve_algo=14”
+
+@#include "purely_backward_common.inc"
+
+perfect_foresight_solver(solve_algo = 14);
+
+ref = load('purely_backward_reference/Output/purely_backward_reference_results.mat');
+
+if max(max(abs(oo_.endo_simul - ref.oo_.endo_simul))) > 5e-5
+ error('Incorrect results for perfect foresight with solve_algo=14')
+end
+
+
diff --git a/tests/solve_algo_12_14/purely_backward_common.inc b/tests/solve_algo_12_14/purely_backward_common.inc
new file mode 100644
index 0000000000000000000000000000000000000000..089bc87df0bf8b0d226dd04629e9bcc678472c8d
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_backward_common.inc
@@ -0,0 +1,34 @@
+@#include "backward_model.inc"
+
+initval;
+@#for i in 1:n
+y@{i} = 1;
+z@{i} = 0.5;
+@#endfor
+end;
+
+shocks;
+@#for i in 1:n
+var e_x@{i};
+periods 1;
+values @{0.1/n};
+
+var e_y@{i};
+periods 2;
+values @{0.05+0.1/n};
+
+var e_z@{i};
+periods 3;
+values @{0.1+0.1/n};
+
+var e_t@{i};
+periods 4;
+values @{0.15+0.1/n};
+
+var e_s@{i};
+periods 5;
+values @{0.2+0.1/n};
+@#endfor
+end;
+
+perfect_foresight_setup(periods = 100);
diff --git a/tests/solve_algo_12_14/purely_backward_reference.mod b/tests/solve_algo_12_14/purely_backward_reference.mod
new file mode 100644
index 0000000000000000000000000000000000000000..397cc4ace1ba4cfd2fca02c38881674aaeb2f3a9
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_backward_reference.mod
@@ -0,0 +1,5 @@
+// Reference perfect foresight simulation of a purely backward model with default “solve_algo” value
+
+@#include "purely_backward_common.inc"
+
+perfect_foresight_solver;
diff --git a/tests/solve_algo_12_14/purely_forward_12.mod b/tests/solve_algo_12_14/purely_forward_12.mod
new file mode 100644
index 0000000000000000000000000000000000000000..46334bf892934ea975912f2dc94e2d50d95ba50f
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_forward_12.mod
@@ -0,0 +1,13 @@
+// Check the correctedness of perfect foresight simulation of a purely forward model with “solve_algo=12”
+
+@#include "purely_forward_common.inc"
+
+perfect_foresight_solver(solve_algo = 12);
+
+ref = load('purely_forward_reference/Output/purely_forward_reference_results.mat');
+
+if max(max(abs(oo_.endo_simul - ref.oo_.endo_simul))) > 2e-4
+ error('Incorrect results for perfect foresight with solve_algo=12')
+end
+
+
diff --git a/tests/solve_algo_12_14/purely_forward_14.mod b/tests/solve_algo_12_14/purely_forward_14.mod
new file mode 100644
index 0000000000000000000000000000000000000000..c3c9eda65bcec8b11d23dd8613a7d8c3b024ecb6
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_forward_14.mod
@@ -0,0 +1,13 @@
+// Check the correctedness of perfect foresight simulation of a purely forward model with “solve_algo=14”
+
+@#include "purely_forward_common.inc"
+
+perfect_foresight_solver(solve_algo = 14);
+
+ref = load('purely_forward_reference/Output/purely_forward_reference_results.mat');
+
+if max(max(abs(oo_.endo_simul - ref.oo_.endo_simul))) > 2e-4
+ error('Incorrect results for perfect foresight with solve_algo=14')
+end
+
+
diff --git a/tests/solve_algo_12_14/purely_forward_common.inc b/tests/solve_algo_12_14/purely_forward_common.inc
new file mode 100644
index 0000000000000000000000000000000000000000..16df0d2dc4c468e05a048bcab7b8b01159f6d8b9
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_forward_common.inc
@@ -0,0 +1,60 @@
+/* A simple purely forward model, consisting of:
+ - a forward VAR(3) with one variable in level, one in log and one log-forward-differentiated;
+ - two other variables that depend on the contemporaneous values of the VAR(3),
+ and which together for a simultaneous block.
+ Those five equations are repeated an arbitrary number of times (n). */
+
+@#define n = 10
+
+@#for i in 1:n
+var x@{i}, y@{i}, z@{i}, t@{i}, s@{i};
+varexo e_x@{i}, e_y@{i}, e_z@{i}, e_t@{i}, e_s@{i};
+@#endfor
+
+model;
+@#for i in 1:n
+x@{i} = 0.1*x@{i}(+1) + 0.2*log(y@{i}(+1)) + 0.3*(log(z@{i})-log(z@{i}(+1))) + e_x@{i};
+[name = 'y']
+log(y@{i}) = 0.4*x@{i}(+1) + 0.5*log(y@{i}(+1)) + 0.6*(log(z@{i})-log(z@{i}(+1))) + e_y@{i};
+[name = 'z']
+/* NB: we choose 0.5 as steady state for z, since initial value is drawn from
+ [0,1] (in the “simul_backward” case) */
+log(z@{i})-log(z@{i}(+1)) = -0.8*(log(z@{i}(+1)) - log(0.5)) + 0.1*x@{i}(+1) + 0.2*log(y@{i}(+1)) - 0.3*(log(z@{i}) - log(z@{i}(+1))) + e_z@{i};
+
+t@{i} = x@{i} + 2*log(y@{i}) + 3*(log(z@{i})-log(z@{i}(+1))) + s@{i} + e_t@{i};
+s@{i} = 4*x@{i}(+1) + 3*t@{i} + e_s@{i};
+@#endfor
+end;
+
+initval;
+@#for i in 1:n
+y@{i} = 1;
+z@{i} = 0.5;
+@#endfor
+end;
+
+shocks;
+@#for i in 1:n
+var e_x@{i};
+periods 99;
+values @{0.1/n};
+
+var e_y@{i};
+periods 98;
+values @{0.05+0.1/n};
+
+var e_z@{i};
+periods 97;
+values @{0.1+0.1/n};
+
+var e_t@{i};
+periods 96;
+values @{0.15+0.1/n};
+
+var e_s@{i};
+periods 95;
+values @{0.2+0.1/n};
+@#endfor
+end;
+
+perfect_foresight_setup(periods = 100);
diff --git a/tests/solve_algo_12_14/purely_forward_reference.mod b/tests/solve_algo_12_14/purely_forward_reference.mod
new file mode 100644
index 0000000000000000000000000000000000000000..70ec7cfb47ba9f06fb8a3582e5e2dd9711319501
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_forward_reference.mod
@@ -0,0 +1,5 @@
+// Reference perfect foresight simulation of a purely forward model with default “solve_algo” value
+
+@#include "purely_forward_common.inc"
+
+perfect_foresight_solver;
diff --git a/tests/solve_algo_12_14/purely_static_12.mod b/tests/solve_algo_12_14/purely_static_12.mod
new file mode 100644
index 0000000000000000000000000000000000000000..9b4da090fea91e4891705a0f03cbf9eb302f3711
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_static_12.mod
@@ -0,0 +1,11 @@
+// Check the correctedness of perfect foresight simulation of a purely static model with “solve_algo=12”
+
+@#include "purely_static_common.inc"
+
+perfect_foresight_solver(solve_algo = 12);
+
+ref = load('purely_static_reference/Output/purely_static_reference_results.mat');
+
+if max(max(abs(oo_.endo_simul - ref.oo_.endo_simul))) > 1e-16
+ error('Incorrect results for perfect foresight with solve_algo=12')
+end
diff --git a/tests/solve_algo_12_14/purely_static_14.mod b/tests/solve_algo_12_14/purely_static_14.mod
new file mode 100644
index 0000000000000000000000000000000000000000..3044f4ebf5e6ab2254ff35dbbc91dbc32bfbda32
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_static_14.mod
@@ -0,0 +1,11 @@
+// Check the correctedness of perfect foresight simulation of a purely static model with “solve_algo=14”
+
+@#include "purely_static_common.inc"
+
+perfect_foresight_solver(solve_algo = 14);
+
+ref = load('purely_static_reference/Output/purely_static_reference_results.mat');
+
+if max(max(abs(oo_.endo_simul - ref.oo_.endo_simul))) > 1e-16
+ error('Incorrect results for perfect foresight with solve_algo=14')
+end
diff --git a/tests/solve_algo_12_14/purely_static_common.inc b/tests/solve_algo_12_14/purely_static_common.inc
new file mode 100644
index 0000000000000000000000000000000000000000..14189814f0ccd9599fa03e6658b4d4bd2808f587
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_static_common.inc
@@ -0,0 +1,32 @@
+@#define n = 10
+
+@#for i in 1:n
+var x@{i}, t@{i}, s@{i};
+varexo e_x@{i}, e_t@{i}, e_s@{i};
+@#endfor
+
+model;
+@#for i in 1:n
+x@{i} = e_x@{i};
+t@{i} = x@{i} + s@{i} + e_t@{i};
+s@{i} = 3*t@{i} + e_s@{i};
+@#endfor
+end;
+
+shocks;
+@#for i in 1:n
+var e_x@{i};
+periods 1;
+values @{0.1/n};
+
+var e_t@{i};
+periods 2;
+values @{0.05+0.1/n};
+
+var e_s@{i};
+periods 3;
+values @{0.1+0.1/n};
+@#endfor
+end;
+
+perfect_foresight_setup(periods = 100);
diff --git a/tests/solve_algo_12_14/purely_static_reference.mod b/tests/solve_algo_12_14/purely_static_reference.mod
new file mode 100644
index 0000000000000000000000000000000000000000..0a619b41c53d245a712188651bf5cd04612b58dc
--- /dev/null
+++ b/tests/solve_algo_12_14/purely_static_reference.mod
@@ -0,0 +1,5 @@
+// Reference perfect foresight simulation of a purely static model with default “solve_algo” value
+
+@#include "purely_static_common.inc"
+
+perfect_foresight_solver;
diff --git a/tests/solve_algo_12_14/simul_backward_12.mod b/tests/solve_algo_12_14/simul_backward_12.mod
new file mode 100644
index 0000000000000000000000000000000000000000..108d8f6039756cf1b5290a3d57e453408687ca2f
--- /dev/null
+++ b/tests/solve_algo_12_14/simul_backward_12.mod
@@ -0,0 +1,12 @@
+// Check the correctedness of “simul_backward” with “solve_algo=12”
+
+@#include "simul_backward_common.inc"
+
+options_.solve_algo = 12;
+s = simul_backward_model(initialconditions, @{simlength}, innovations);
+
+ref = load('simul_backward_ref.mat');
+
+if max(max(abs(s.data - ref.s.data))) > 5e-4
+ error('Incorrect results for simul_backward with solve_algo=12')
+end
diff --git a/tests/solve_algo_12_14/simul_backward_14.mod b/tests/solve_algo_12_14/simul_backward_14.mod
new file mode 100644
index 0000000000000000000000000000000000000000..8816f936032c4a1578f764689a13299a7d6cd9ef
--- /dev/null
+++ b/tests/solve_algo_12_14/simul_backward_14.mod
@@ -0,0 +1,12 @@
+// Check the correctedness of “simul_backward” with “solve_algo=14”
+
+@#include "simul_backward_common.inc"
+
+options_.solve_algo = 14;
+s = simul_backward_model(initialconditions, @{simlength}, innovations);
+
+ref = load('simul_backward_ref.mat');
+
+if max(max(abs(s.data - ref.s.data))) > 5e-4
+ error('Incorrect results for simul_backward with solve_algo=14')
+end
diff --git a/tests/solve_algo_12_14/simul_backward_common.inc b/tests/solve_algo_12_14/simul_backward_common.inc
new file mode 100644
index 0000000000000000000000000000000000000000..86f7fabccba54365f74c7e3fa9a17463637eceb7
--- /dev/null
+++ b/tests/solve_algo_12_14/simul_backward_common.inc
@@ -0,0 +1,15 @@
+@#include "backward_model.inc"
+
+@#define simlength = 100
+
+initialconditions = dseries(rand(2, @{5*n}), '2021Q3', { ...
+@#for i in 1:n
+ 'x@{i}', 'y@{i}', 'z@{i}', 't@{i}', 's@{i}', ...
+@#endfor
+});
+
+innovations = dseries(randn(@{simlength}, @{5*n}), '2022Q1', { ...
+@#for i in 1:n
+ 'e_x@{i}', 'e_y@{i}', 'e_z@{i}', 'e_t@{i}', 'e_s@{i}', ...
+@#endfor
+});
diff --git a/tests/solve_algo_12_14/simul_backward_reference.mod b/tests/solve_algo_12_14/simul_backward_reference.mod
new file mode 100644
index 0000000000000000000000000000000000000000..cf5284140d87b9fcf6c6f7fa51b0f8a631845485
--- /dev/null
+++ b/tests/solve_algo_12_14/simul_backward_reference.mod
@@ -0,0 +1,7 @@
+// Reference simulation with “simul_backward” and default “solve_algo” value
+
+@#include "simul_backward_common.inc"
+
+s = simul_backward_model(initialconditions, @{simlength}, innovations);
+
+save simul_backward_ref.mat s