Design and performance improvement to solve_algo={12,14}
Use the new time-recursive block decomposition computed by the preprocessor for: - the simulation of backward models with “simul_backward” - the perfect foresight simulation of purely backward/forward/static models Also note that in this case, the preprocessor now defaults to “mfs=3” (i.e. it minimizes the set of feedback variables and tries to renormalize equations). This replaces the previous algorithm based on Dulmage-Mendelsohn (dmperm), plus an ad hoc identification of some equations that can be evaluated (those with a LHS equal to a variable, the log of a variable, or the diff-log of a variable). By the way, the block_trust_region MEX has been modified so that it accepts a boolean argument to decide whether it performs a Dulmage-Mendelsohn decomposition (if not, then it performs a simple trust region on the whole nonlinear system). This provides a significant performance improvement (of almost an order of magnitude for solve_algo=14 on a 700 equations model).unit-tests
Sébastien Villemot
parent
10fdc42516
commit
d574705b4a
|
|
@ -3011,27 +3011,27 @@ Finding the steady state with Dynare nonlinear solver
|
||||||
|
|
||||||
``12``
|
``12``
|
||||||
|
|
||||||
Specialized version of ``2`` for models where all the equations
|
Computes a block decomposition and then applies a Newton-type
|
||||||
have one endogenous variable on the left hand side and where
|
solver on those smaller blocks rather than on the full
|
||||||
each equation determines a different endogenous variable. Only
|
nonlinear system. This is similar to ``2``, but is typically
|
||||||
expressions allowed on the left hand side are the natural
|
more efficient. The block decomposition is done at the
|
||||||
logarithm of an endogenous variable, the first difference of an
|
preprocessor level, which brings two benefits: it identifies
|
||||||
endogenous variable (with the ``diff`` operator), or the first
|
blocks that can be evaluated rather than solved; and evulations
|
||||||
difference of the logarithm of an endogenous variable.
|
of the residual and Jacobian of the model are more efficient
|
||||||
Univariate blocks are solved by evaluating the expression on the
|
because only the relevant elements are recomputed at every
|
||||||
right hand side.
|
iteration.
|
||||||
|
|
||||||
``14``
|
``14``
|
||||||
|
|
||||||
Specialized version of ``4`` for models where all the equations
|
Computes a block decomposition and then applies a trust region
|
||||||
have one endogenous variable on the left hand side and where
|
solver with autoscaling on those smaller blocks rather than on
|
||||||
each equation determines a different endogenous variable. Only
|
the full nonlinear system. This is similar to ``4``, but is
|
||||||
expressions allowed on the left hand side are the natural
|
typically more efficient. The block decomposition is done at
|
||||||
logarithm of an endogenous variable, the first difference of an
|
the preprocessor level, which brings two benefits: it
|
||||||
endogenous variable (with the ``diff`` operator), or the first
|
identifies blocks that can be evaluated rather than solved; and
|
||||||
difference of the logarithm of an endogenous variable..
|
evulations of the residual and Jacobian of the model are more
|
||||||
Univariate blocks are solved by evaluating the expression on the
|
efficient because only the relevant elements are recomputed at
|
||||||
right hand side.
|
every iteration.
|
||||||
|
|
||||||
|br| Default value is ``4``.
|
|br| Default value is ``4``.
|
||||||
|
|
||||||
|
|
|
||||||
|
|
@ -47,10 +47,6 @@ if ~M_.maximum_lag
|
||||||
return
|
return
|
||||||
end
|
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
|
if nargin<3
|
||||||
Innovations = [];
|
Innovations = [];
|
||||||
else
|
else
|
||||||
|
|
|
||||||
|
|
@ -40,36 +40,64 @@ function [ysim, xsim] = simul_backward_nonlinear_model_(initialconditions, sampl
|
||||||
|
|
||||||
debug = false;
|
debug = false;
|
||||||
|
|
||||||
model_dynamic_s = str2func('dynamic_backward_model_for_simulation');
|
|
||||||
|
|
||||||
if ~isempty(innovations)
|
if ~isempty(innovations)
|
||||||
DynareOutput.exo_simul(initialconditions.nobs+(1:samplesize),:) = innovations;
|
DynareOutput.exo_simul(initialconditions.nobs+(1:samplesize),:) = innovations;
|
||||||
end
|
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).
|
% Simulations (call a Newton-like algorithm for each period).
|
||||||
for it = initialconditions.nobs+(1:samplesize)
|
for it = initialconditions.nobs+(1:samplesize)
|
||||||
if debug
|
if debug
|
||||||
dprintf('Period t = %s.', num2str(it-initialconditions.nobs));
|
dprintf('Period t = %s.', num2str(it-initialconditions.nobs));
|
||||||
end
|
end
|
||||||
y_ = DynareOutput.endo_simul(:,it-1);
|
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.
|
y = y_; % A good guess for the initial conditions is the previous values for the endogenous variables.
|
||||||
try
|
try
|
||||||
if ismember(DynareOptions.solve_algo, [12,14])
|
if ismember(DynareOptions.solve_algo, [12,14])
|
||||||
[DynareOutput.endo_simul(:,it), errorflag, ~, ~, errorcode] = dynare_solve(model_dynamic_s, y, ...
|
x = DynareOutput.exo_simul(it,:);
|
||||||
DynareOptions.simul.maxit, DynareOptions.dynatol.f, DynareOptions.dynatol.x, ...
|
T = NaN(DynareModel.block_structure.dyn_tmp_nbr);
|
||||||
DynareOptions, ...
|
y_dynamic = [y_; y; NaN(DynareModel.endo_nbr, 1)];
|
||||||
DynareModel.isloggedlhs, DynareModel.isauxdiffloggedrhs, DynareModel.endo_names, DynareModel.lhs, ...
|
for blk = 1:length(DynareModel.block_structure.block)
|
||||||
model_dynamic, ylag, DynareOutput.exo_simul, DynareModel.params, DynareOutput.steady_state, it);
|
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
|
else
|
||||||
[DynareOutput.endo_simul(:,it), errorflag, ~, ~, errorcode] = ...
|
[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.simul.maxit, DynareOptions.dynatol.f, DynareOptions.dynatol.x, ...
|
||||||
DynareOptions, ...
|
DynareOptions, ...
|
||||||
model_dynamic, ylag, DynareOutput.exo_simul, DynareModel.params, DynareOutput.steady_state, it);
|
model_dynamic, y_(iy1), DynareOutput.exo_simul, DynareModel.params, DynareOutput.steady_state, it);
|
||||||
end
|
if errorflag
|
||||||
if errorflag
|
error('Nonlinear solver routine failed with errorcode=%i in period %i.', errorcode, it)
|
||||||
error('Nonlinear solver routine failed with errorcode=%i in period %i.', errorcode, it)
|
end
|
||||||
end
|
end
|
||||||
catch
|
catch
|
||||||
DynareOutput.endo_simul = DynareOutput.endo_simul(:, 1:it-1);
|
DynareOutput.endo_simul = DynareOutput.endo_simul(:, 1:it-1);
|
||||||
|
|
@ -111,7 +139,7 @@ for it = initialconditions.nobs+(1:samplesize)
|
||||||
%
|
%
|
||||||
% Evaluate and check the residuals
|
% 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_nan = isnan(r);
|
||||||
residuals_evaluating_to_inf = isinf(r);
|
residuals_evaluating_to_inf = isinf(r);
|
||||||
residuals_evaluating_to_complex = ~isreal(r);
|
residuals_evaluating_to_complex = ~isreal(r);
|
||||||
|
|
@ -142,6 +170,8 @@ end
|
||||||
ysim = DynareOutput.endo_simul(1:DynareModel.orig_endo_nbr,:);
|
ysim = DynareOutput.endo_simul(1:DynareModel.orig_endo_nbr,:);
|
||||||
xsim = DynareOutput.exo_simul;
|
xsim = DynareOutput.exo_simul;
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
function display_names_of_problematic_equations(DynareModel, eqtags, TruthTable)
|
function display_names_of_problematic_equations(DynareModel, eqtags, TruthTable)
|
||||||
for i=1:DynareModel.orig_endo_nbr
|
for i=1:DynareModel.orig_endo_nbr
|
||||||
if TruthTable(i)
|
if TruthTable(i)
|
||||||
|
|
@ -153,3 +183,4 @@ for i=DynareModel.orig_endo_nbr+1:DynareModel.endo_nbr
|
||||||
dprintf(' - Auxiliary equation for %s', DynareModel.endo_names{i})
|
dprintf(' - Auxiliary equation for %s', DynareModel.endo_names{i})
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
end
|
||||||
|
|
|
||||||
|
|
@ -58,21 +58,10 @@ else
|
||||||
in0 = nn;
|
in0 = nn;
|
||||||
end
|
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
|
% checking initial values
|
||||||
% TODO We should have an option to deactivate the randomization.
|
% TODO We should have an option to deactivate the randomization.
|
||||||
if jacobian_flag
|
if jacobian_flag
|
||||||
[fvec, fjac] = feval(f, x, args{:});
|
[fvec, fjac] = feval(f, x, varargin{:});
|
||||||
wrong_initial_guess_flag = false;
|
wrong_initial_guess_flag = false;
|
||||||
if ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:)))) || any(~isreal(fvec)) || any(~isreal(fjac(:)))
|
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
|
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
|
while wrong_initial_guess_flag && tentative_number<=in0*10
|
||||||
tentative_number = tentative_number+1;
|
tentative_number = tentative_number+1;
|
||||||
x(idx) = rand(in0, 1)*10;
|
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(:))));
|
wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:))));
|
||||||
end
|
end
|
||||||
% If all previous attempts failed, try with real numbers.
|
% 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
|
while wrong_initial_guess_flag && tentative_number<=in0*10
|
||||||
tentative_number = tentative_number+1;
|
tentative_number = tentative_number+1;
|
||||||
x(idx) = randn(in0, 1)*10;
|
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(:))));
|
wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:))));
|
||||||
end
|
end
|
||||||
% Last tentative, ff all previous attempts failed, try with negative numbers.
|
% 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
|
while wrong_initial_guess_flag && tentative_number<=in0*10
|
||||||
tentative_number = tentative_number+1;
|
tentative_number = tentative_number+1;
|
||||||
x(idx) = -rand(in0, 1)*10;
|
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(:))));
|
wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:))));
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
else
|
else
|
||||||
fvec = feval(f, x, args{:});
|
fvec = feval(f, x, varargin{:});
|
||||||
fjac = zeros(nn, nn);
|
fjac = zeros(nn, nn);
|
||||||
if ~ismember(options.solve_algo,[10,11]) && max(abs(fvec)) < tolf
|
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)
|
% 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
|
while wrong_initial_guess_flag && tentative_number<=in0*10
|
||||||
tentative_number = tentative_number+1;
|
tentative_number = tentative_number+1;
|
||||||
x(idx) = rand(in0, 1)*10;
|
x(idx) = rand(in0, 1)*10;
|
||||||
fvec = feval(f, x, args{:});
|
fvec = feval(f, x, varargin{:});
|
||||||
wrong_initial_guess_flag = ~all(isfinite(fvec));
|
wrong_initial_guess_flag = ~all(isfinite(fvec));
|
||||||
end
|
end
|
||||||
% If all previous attempts failed, try with real numbers.
|
% If all previous attempts failed, try with real numbers.
|
||||||
|
|
@ -133,7 +122,7 @@ else
|
||||||
while wrong_initial_guess_flag && tentative_number<=in0*10
|
while wrong_initial_guess_flag && tentative_number<=in0*10
|
||||||
tentative_number = tentative_number+1;
|
tentative_number = tentative_number+1;
|
||||||
x(idx) = randn(in0, 1)*10;
|
x(idx) = randn(in0, 1)*10;
|
||||||
fvec = feval(f, x, args{:});
|
fvec = feval(f, x, varargin{:});
|
||||||
wrong_initial_guess_flag = ~all(isfinite(fvec));
|
wrong_initial_guess_flag = ~all(isfinite(fvec));
|
||||||
end
|
end
|
||||||
% Last tentative, ff all previous attempts failed, try with negative numbers.
|
% 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
|
while wrong_initial_guess_flag && tentative_number<=in0*10
|
||||||
tentative_number = tentative_number+1;
|
tentative_number = tentative_number+1;
|
||||||
x(idx) = -rand(in0, 1)*10;
|
x(idx) = -rand(in0, 1)*10;
|
||||||
fvec = feval(f, x, args{:});
|
fvec = feval(f, x, varargin{:});
|
||||||
wrong_initial_guess_flag = ~all(isfinite(fvec));
|
wrong_initial_guess_flag = ~all(isfinite(fvec));
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
@ -199,7 +188,7 @@ if options.solve_algo == 0
|
||||||
end
|
end
|
||||||
|
|
||||||
if ~isoctave
|
if ~isoctave
|
||||||
[x, fvec, errorcode, ~, fjac] = fsolve(f, x, options4fsolve, args{:});
|
[x, fvec, errorcode, ~, fjac] = fsolve(f, x, options4fsolve, varargin{:});
|
||||||
else
|
else
|
||||||
% Under Octave, use a wrapper, since fsolve() does not have a 4th arg
|
% Under Octave, use a wrapper, since fsolve() does not have a 4th arg
|
||||||
if ischar(f)
|
if ischar(f)
|
||||||
|
|
@ -207,7 +196,7 @@ if options.solve_algo == 0
|
||||||
else
|
else
|
||||||
f2 = f;
|
f2 = f;
|
||||||
end
|
end
|
||||||
[x, fvec, errorcode, ~, fjac] = fsolve(@(x) f2(x, args{:}), x, options4fsolve);
|
[x, fvec, errorcode, ~, fjac] = fsolve(@(x) f2(x, varargin{:}), x, options4fsolve);
|
||||||
end
|
end
|
||||||
if errorcode==1
|
if errorcode==1
|
||||||
errorflag = false;
|
errorflag = false;
|
||||||
|
|
@ -220,91 +209,42 @@ if options.solve_algo == 0
|
||||||
else
|
else
|
||||||
errorflag = true;
|
errorflag = true;
|
||||||
end
|
end
|
||||||
elseif options.solve_algo==1
|
elseif ismember(options.solve_algo, [1, 12])
|
||||||
[x, errorflag, errorcode] = solve1(f, x, 1:nn, 1:nn, jacobian_flag, options.gstep, tolf, tolx, maxit, [], options.debug, args{:});
|
%% NB: It is the responsibility of the caller to deal with the block decomposition if solve_algo=12
|
||||||
[fvec, fjac] = feval(f, x, args{:});
|
[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
|
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{:});
|
[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, args{:});
|
[fvec, fjac] = feval(f, x, varargin{:});
|
||||||
elseif ismember(options.solve_algo, [2, 12, 4])
|
elseif ismember(options.solve_algo, [2, 4])
|
||||||
if ismember(options.solve_algo, [2, 12])
|
if options.solve_algo == 2
|
||||||
solver = @solve1;
|
solver = @solve1;
|
||||||
else
|
else
|
||||||
solver = @trust_region;
|
solver = @trust_region;
|
||||||
end
|
end
|
||||||
specializedunivariateblocks = options.solve_algo == 12;
|
|
||||||
if ~jacobian_flag
|
if ~jacobian_flag
|
||||||
fjac = zeros(nn,nn) ;
|
fjac = zeros(nn,nn) ;
|
||||||
dh = max(abs(x), options.gstep(1)*ones(nn,1))*eps^(1/3);
|
dh = max(abs(x), options.gstep(1)*ones(nn,1))*eps^(1/3);
|
||||||
for j = 1:nn
|
for j = 1:nn
|
||||||
xdh = x ;
|
xdh = x ;
|
||||||
xdh(j) = xdh(j)+dh(j) ;
|
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
|
||||||
end
|
end
|
||||||
[j1,j2,r,s] = dmperm(fjac);
|
[j1,j2,r,s] = dmperm(fjac);
|
||||||
JAC = abs(fjac(j1,j2))>0;
|
|
||||||
if options.debug
|
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
|
end
|
||||||
l = 0;
|
|
||||||
fre = false;
|
|
||||||
for i=length(r)-1:-1:1
|
for i=length(r)-1:-1:1
|
||||||
blocklength = r(i+1)-r(i);
|
blocklength = r(i+1)-r(i);
|
||||||
if options.debug
|
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
|
end
|
||||||
j = r(i):r(i+1)-1;
|
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);
|
blockcolumns=s(i+1)-s(i);
|
||||||
if blockcolumns ~= blocklength
|
if blockcolumns ~= blocklength
|
||||||
%non-square-block in DM; check whether initial value is solution
|
%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
|
if norm(fval_check(j1(j))) < tolf
|
||||||
errorflag = false;
|
errorflag = false;
|
||||||
errorcode = 0;
|
errorcode = 0;
|
||||||
|
|
@ -317,10 +257,9 @@ elseif ismember(options.solve_algo, [2, 12, 4])
|
||||||
options.gstep, ...
|
options.gstep, ...
|
||||||
tolf, options.solve_tolx, maxit, ...
|
tolf, options.solve_tolx, maxit, ...
|
||||||
options.trust_region_initial_step_bound_factor, ...
|
options.trust_region_initial_step_bound_factor, ...
|
||||||
options.debug, args{:});
|
options.debug, varargin{:});
|
||||||
fre = true;
|
|
||||||
else
|
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
|
%overdetermined block
|
||||||
errorflag = true;
|
errorflag = true;
|
||||||
errorcode = 0;
|
errorcode = 0;
|
||||||
|
|
@ -329,31 +268,31 @@ elseif ismember(options.solve_algo, [2, 12, 4])
|
||||||
return
|
return
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
fvec = feval(f, x, args{:});
|
fvec = feval(f, x, varargin{:});
|
||||||
if max(abs(fvec))>tolf
|
if max(abs(fvec))>tolf
|
||||||
disp_verbose('Call solver on the full nonlinear problem.',options.verbosity)
|
disp_verbose('Call solver on the full nonlinear problem.',options.verbosity)
|
||||||
[x, errorflag, errorcode] = solver(f, x, 1:nn, 1:nn, jacobian_flag, ...
|
[x, errorflag, errorcode] = solver(f, x, 1:nn, 1:nn, jacobian_flag, ...
|
||||||
options.gstep, tolf, options.solve_tolx, maxit, ...
|
options.gstep, tolf, options.solve_tolx, maxit, ...
|
||||||
options.trust_region_initial_step_bound_factor, ...
|
options.trust_region_initial_step_bound_factor, ...
|
||||||
options.debug, args{:});
|
options.debug, varargin{:});
|
||||||
end
|
end
|
||||||
[fvec, fjac] = feval(f, x, args{:});
|
[fvec, fjac] = feval(f, x, varargin{:});
|
||||||
elseif options.solve_algo==3
|
elseif options.solve_algo==3
|
||||||
if jacobian_flag
|
if jacobian_flag
|
||||||
[x, errorcode] = csolve(f, x, f, tolf, maxit, args{:});
|
[x, errorcode] = csolve(f, x, f, tolf, maxit, varargin{:});
|
||||||
else
|
else
|
||||||
[x, errorcode] = csolve(f, x, [], tolf, maxit, args{:});
|
[x, errorcode] = csolve(f, x, [], tolf, maxit, varargin{:});
|
||||||
end
|
end
|
||||||
if errorcode==0
|
if errorcode==0
|
||||||
errorflag = false;
|
errorflag = false;
|
||||||
else
|
else
|
||||||
errorflag = true;
|
errorflag = true;
|
||||||
end
|
end
|
||||||
[fvec, fjac] = feval(f, x, args{:});
|
[fvec, fjac] = feval(f, x, varargin{:});
|
||||||
elseif options.solve_algo==10
|
elseif options.solve_algo==10
|
||||||
% LMMCP
|
% LMMCP
|
||||||
olmmcp = options.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_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);
|
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;
|
fvec(eq_to_ignore)=0;
|
||||||
|
|
@ -373,7 +312,7 @@ elseif options.solve_algo == 11
|
||||||
omcppath = options.mcppath;
|
omcppath = options.mcppath;
|
||||||
global mcp_data
|
global mcp_data
|
||||||
mcp_data.func = f;
|
mcp_data.func = f;
|
||||||
mcp_data.args = args;
|
mcp_data.args = varargin;
|
||||||
try
|
try
|
||||||
[x, fval, jac, mu] = pathmcp(x,omcppath.lb,omcppath.ub,'mcp_func',omcppath.A,omcppath.b,omcppath.t,omcppath.mu0);
|
[x, fval, jac, mu] = pathmcp(x,omcppath.lb,omcppath.ub,'mcp_func',omcppath.A,omcppath.b,omcppath.t,omcppath.mu0);
|
||||||
catch
|
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);
|
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;
|
fvec(eq_to_ignore)=0;
|
||||||
elseif ismember(options.solve_algo, [13, 14])
|
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
|
if ~jacobian_flag
|
||||||
error('DYNARE_SOLVE: option solve_algo=13|14 needs computed Jacobian')
|
error('DYNARE_SOLVE: option solve_algo=13 needs computed Jacobian')
|
||||||
end
|
end
|
||||||
auxstruct = struct();
|
[x, errorflag, errorcode] = block_trust_region(f, x, tolf, options.solve_tolx, maxit, ...
|
||||||
if options.solve_algo == 14
|
options.trust_region_initial_step_bound_factor, ...
|
||||||
auxstruct.lhs = lhs;
|
options.solve_algo == 13, ... % Only block-decompose with Dulmage-Mendelsohn for 13, not for 14
|
||||||
auxstruct.endo_names = endo_names;
|
options.debug, varargin{:});
|
||||||
auxstruct.isloggedlhs = isloggedlhs;
|
[fvec, fjac] = feval(f, x, varargin{:});
|
||||||
auxstruct.isauxdiffloggedrhs = isauxdiffloggedrhs;
|
|
||||||
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{:});
|
|
||||||
else
|
else
|
||||||
error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,9,10,11,12,13,14]')
|
error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,9,10,11,12,13,14]')
|
||||||
end
|
end
|
||||||
|
|
|
||||||
|
|
@ -50,6 +50,9 @@ if options_.linear && (isequal(options_.stack_solve_algo, 0) || isequal(options_
|
||||||
end
|
end
|
||||||
|
|
||||||
if options_.block
|
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
|
if options_.bytecode
|
||||||
try
|
try
|
||||||
oo_.endo_simul = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state,1, periods+2), periods);
|
oo_.endo_simul = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state,1, periods+2), periods);
|
||||||
|
|
|
||||||
|
|
@ -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
|
||||||
|
|
@ -20,38 +20,71 @@ function [endogenousvariables, info] = sim1_purely_backward(endogenousvariables,
|
||||||
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
ny0 = nnz(M.lead_lag_incidence(2,:)); % Number of variables at current period
|
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
|
if ny0 ~= M.endo_nbr
|
||||||
error('All endogenous variables must appear at the current period!')
|
error('All endogenous variables must appear at the current period!')
|
||||||
end
|
end
|
||||||
|
|
||||||
if ismember(options.solve_algo, [12,14]) && ~M.possible_to_use_solve_algo_12_14
|
if ismember(options.solve_algo, [12,14])
|
||||||
error(M.message_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
|
end
|
||||||
|
|
||||||
dynamicmodel = str2func(sprintf('%s.%s', M.fname, 'dynamic'));
|
function [r, J] = block_wrapper(z, feedback_vars_idx, func, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T)
|
||||||
dynamicmodel_s = str2func('dynamic_backward_model_for_simulation');
|
% 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;
|
info.status = true;
|
||||||
|
|
||||||
for it = M.maximum_lag + (1:options.periods)
|
for it = M.maximum_lag + (1:options.periods)
|
||||||
y = endogenousvariables(:,it-1); % Values at previous period, also used as guess value for current period
|
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])
|
if ismember(options.solve_algo, [12,14])
|
||||||
[tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, y, ...
|
x = exogenousvariables(it,:);
|
||||||
options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
|
T = NaN(M.block_structure.dyn_tmp_nbr);
|
||||||
options, M.isloggedlhs, M.isauxdiffloggedrhs, M.endo_names, M.lhs, ...
|
y_dynamic = [y; y; NaN(M.endo_nbr, 1)];
|
||||||
dynamicmodel, ylag, exogenousvariables, M.params, steadystate, it);
|
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
|
else
|
||||||
|
ylag = y(iyb);
|
||||||
[tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, y, ...
|
[tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, y, ...
|
||||||
options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
|
options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
|
||||||
options, dynamicmodel, ylag, exogenousvariables, M.params, steadystate, it);
|
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
|
end
|
||||||
if check
|
end
|
||||||
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
|
end
|
||||||
|
|
@ -19,34 +19,71 @@ function [endogenousvariables, info] = sim1_purely_forward(endogenousvariables,
|
||||||
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
ny0 = nnz(M.lead_lag_incidence(1,:)); % Number of variables at current period
|
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
|
if ny0 ~= M.endo_nbr
|
||||||
error('All endogenous variables must appear at the current period!')
|
error('All endogenous variables must appear at the current period!')
|
||||||
end
|
end
|
||||||
|
|
||||||
dynamicmodel = str2func(sprintf('%s.%s', M.fname, 'dynamic'));
|
if ismember(options.solve_algo, [12,14])
|
||||||
dynamicmodel_s = str2func('dynamic_forward_model_for_simulation');
|
[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;
|
info.status = true;
|
||||||
|
|
||||||
for it = options.periods:-1:1
|
for it = options.periods:-1:1
|
||||||
yf = endogenousvariables(:,it+1); % Values at next period, also used as guess value for current period
|
yf = endogenousvariables(:,it+1); % Values at next period, also used as guess value for current period
|
||||||
if ismember(options.solve_algo, [12,14])
|
if ismember(options.solve_algo, [12,14])
|
||||||
[tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, yf, ...
|
x = exogenousvariables(it,:);
|
||||||
options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
|
T = NaN(M.block_structure.dyn_tmp_nbr);
|
||||||
options, M.isloggedlhs, M.isauxdiffloggedrhs, M.endo_names, M.lhs, ...
|
y_dynamic = [NaN(M.endo_nbr, 1); yf; yf];
|
||||||
dynamicmodel, yf(iyf), exogenousvariables, M.params, steadystate, it);
|
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
|
else
|
||||||
[tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, yf, ...
|
[tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, yf, ...
|
||||||
options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
|
options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
|
||||||
options, ...
|
options, ...
|
||||||
dynamicmodel, yf(iyf), exogenousvariables, M.params, steadystate, it);
|
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
|
end
|
||||||
if check
|
end
|
||||||
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
|
end
|
||||||
|
|
@ -23,12 +23,20 @@ if nnz(M.lead_lag_incidence(1,:)) ~= M.endo_nbr
|
||||||
error('All endogenous variables must appear at the current period!')
|
error('All endogenous variables must appear at the current period!')
|
||||||
end
|
end
|
||||||
|
|
||||||
if ismember(options.solve_algo, [12,14]) && ~M.possible_to_use_solve_algo_12_14
|
if ismember(options.solve_algo, [12,14])
|
||||||
error(M.message_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
|
end
|
||||||
|
|
||||||
dynamicmodel = str2func(sprintf('%s.%s', M.fname, 'dynamic'));
|
function [r, J] = block_wrapper(z, feedback_vars_idx, func, y_dynamic, x, sparse_rowval, sparse_colval, sparse_colptr, T)
|
||||||
dynamicmodel_s = str2func('dynamic_static_model_for_simulation');
|
% 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;
|
info.status = true;
|
||||||
|
|
||||||
|
|
@ -36,21 +44,46 @@ y = endogenousvariables(:,1);
|
||||||
|
|
||||||
for it = 1:options.periods
|
for it = 1:options.periods
|
||||||
if ismember(options.solve_algo, [12,14])
|
if ismember(options.solve_algo, [12,14])
|
||||||
[tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, y, ...
|
x = exogenousvariables(it,:);
|
||||||
options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
|
T = NaN(M.block_structure.dyn_tmp_nbr);
|
||||||
options, M.isloggedlhs, M.isauxdiffloggedrhs, M.endo_names, M.lhs, ...
|
y_dynamic = [NaN(M.endo_nbr, 1); y; NaN(M.endo_nbr, 1)];
|
||||||
dynamicmodel, exogenousvariables, M.params, steadystate, it);
|
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
|
else
|
||||||
[tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, y, ...
|
[tmp, check, ~, ~, errorcode] = dynare_solve(dynamicmodel_s, y, ...
|
||||||
options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
|
options.simul.maxit, options.dynatol.f, options.dynatol.x, ...
|
||||||
options, dynamicmodel, exogenousvariables, M.params, steadystate, it);
|
options, dynamicmodel, exogenousvariables, M.params, steadystate, it);
|
||||||
end
|
if check
|
||||||
if check
|
info.status = false;
|
||||||
info.status = false;
|
if options.debug
|
||||||
if options.debug
|
dprintf('sim1_purely_static: Nonlinear solver routine failed with errorcode=%i in period %i.', errorcode, it)
|
||||||
dprintf('sim1_purely_static: Nonlinear solver routine failed with errorcode=%i in period %i.', errorcode, it)
|
end
|
||||||
end
|
end
|
||||||
|
endogenousvariables(:,it) = tmp;
|
||||||
end
|
end
|
||||||
endogenousvariables(:,it) = tmp;
|
|
||||||
y = endogenousvariables(:,it);
|
y = endogenousvariables(:,it);
|
||||||
end
|
end
|
||||||
|
|
||||||
|
end
|
||||||
|
|
|
||||||
|
|
@ -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
|
|
||||||
|
|
@ -35,16 +35,11 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
|
||||||
integer :: maxiter
|
integer :: maxiter
|
||||||
real(real64), dimension(:), allocatable :: fvec
|
real(real64), dimension(:), allocatable :: fvec
|
||||||
real(real64), dimension(:,:), allocatable :: fjac
|
real(real64), dimension(:,:), allocatable :: fjac
|
||||||
logical :: debug, specializedunivariateblocks
|
logical :: debug, block_decompose
|
||||||
character(len=80) :: debug_msg
|
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
|
if (nrhs < 8 .or. nlhs /= 3) then
|
||||||
call mexErrMsgTxt("Must have at least 7 inputs and exactly 3 outputs")
|
call mexErrMsgTxt("Must have at least 8 inputs and exactly 3 outputs")
|
||||||
end if
|
end if
|
||||||
|
|
||||||
if (.not. ((mxIsChar(prhs(1)) .and. mxGetM(prhs(1)) == 1) .or. mxIsClass(prhs(1), "function_handle"))) then
|
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
|
end if
|
||||||
|
|
||||||
if (.not. (mxIsLogicalScalar(prhs(7)))) then
|
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
|
end if
|
||||||
|
|
||||||
if (.not. (mxIsStruct(prhs(8)) .and. &
|
if (.not. (mxIsLogicalScalar(prhs(8)))) then
|
||||||
(mxGetNumberOfFields(prhs(8)) == 0 .or. mxGetNumberOfFields(prhs(8)) == 4))) then
|
call mexErrMsgTxt("Eigth argument (debug) should be a logical scalar")
|
||||||
call mexErrMsgTxt("Eighth argument should be a struct with either 0 or 4 fields")
|
|
||||||
end if
|
end if
|
||||||
|
|
||||||
specializedunivariateblocks = (mxGetNumberOfFields(prhs(8)) == 4)
|
|
||||||
|
|
||||||
func => prhs(1)
|
func => prhs(1)
|
||||||
tolf = mxGetScalar(prhs(3))
|
tolf = mxGetScalar(prhs(3))
|
||||||
tolx = mxGetScalar(prhs(4))
|
tolx = mxGetScalar(prhs(4))
|
||||||
maxiter = int(mxGetScalar(prhs(5)))
|
maxiter = int(mxGetScalar(prhs(5)))
|
||||||
factor = mxGetScalar(prhs(6))
|
factor = mxGetScalar(prhs(6))
|
||||||
debug = mxGetScalar(prhs(7)) == 1._c_double
|
block_decompose = mxGetScalar(prhs(7)) == 1._c_double
|
||||||
extra_args => prhs(9:nrhs) ! Extra arguments to func are in argument 9 and subsequent ones
|
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)))
|
associate (x_mat => mxGetPr(prhs(2)))
|
||||||
allocate(x(size(x_mat)))
|
allocate(x(size(x_mat)))
|
||||||
x = x_mat
|
x = x_mat
|
||||||
end associate
|
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(fvec(size(x)))
|
||||||
allocate(fjac(size(x), size(x)))
|
allocate(fjac(size(x), size(x)))
|
||||||
|
|
||||||
! Compute block decomposition
|
if (block_decompose) then
|
||||||
nullify(x_indices, f_indices, x_all)
|
! Compute block decomposition
|
||||||
call matlab_fcn(x, fvec, fjac)
|
nullify(x_indices, f_indices, x_all)
|
||||||
call dm_blocks(fjac, blocks)
|
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
|
|
||||||
|
|
||||||
! Solve the system, starting from bottom-rightmost block
|
|
||||||
do i = size(blocks),1,-1
|
|
||||||
if (debug) then
|
if (debug) then
|
||||||
write (debug_msg, "('DYNARE_SOLVE (solve_algo=13|14): solving block ', i0, ' of size ', i0)") &
|
write (debug_msg, "('DYNARE_SOLVE (solve_algo=13): number of blocks = ', i0)") size(blocks)
|
||||||
i, size(blocks(i)%col_indices)
|
|
||||||
call mexPrintf_trim_newline(debug_msg)
|
call mexPrintf_trim_newline(debug_msg)
|
||||||
end if
|
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
|
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)
|
call mexPrintf_trim_newline(debug_msg)
|
||||||
end if
|
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
|
if (size(blocks(i)%col_indices) /= size(blocks(i)%row_indices)) then
|
||||||
! Reevaluation of the residuals is required because the current RHS depends on
|
! Non-square block in DM decomposition
|
||||||
! variables that potentially have been updated previously.
|
! Before erroring out, check whether we are not already at the solution for this block
|
||||||
nullify(x_indices, f_indices, x_all)
|
! See also #1851
|
||||||
call matlab_fcn(x, fvec)
|
if (norm2(fvec(blocks(i)%row_indices)) < tolf) then
|
||||||
deallocate(evaled_cols) ! This shouldn’t be necessary, but it crashes otherwise with gfortran 8
|
cycle
|
||||||
allocate(evaled_cols(0))
|
else
|
||||||
fre = .false.
|
call mexErrMsgTxt("DYNARE_SOLVE (solve_algo=13): the Dulmage-Mendelsohn &
|
||||||
end if
|
&decomposition returned a non-square block. This means that the &
|
||||||
evaled_cols = [ evaled_cols, var]
|
&Jacobian is singular. You may want to try another value for solve_algo.")
|
||||||
block
|
end if
|
||||||
! An associate() construct for lhs_eq and endo_name_var makes the
|
end if
|
||||||
! code crash (with double free) using gfortran 8. Hence use a block
|
|
||||||
character(kind=c_char, len=:), allocatable :: lhs_eq, endo_name_var
|
block
|
||||||
lhs_eq = mxArrayToString(mxGetCell(lhs, int(eq, mwIndex)))
|
real(real64), dimension(size(blocks(i)%col_indices)) :: x_block
|
||||||
endo_name_var = mxArrayToString(mxGetCell(endo_names, int(var, mwIndex)))
|
x_indices => blocks(i)%col_indices
|
||||||
if (lhs_eq == endo_name_var .or. lhs_eq == "log(" // endo_name_var // ")") then
|
f_indices => blocks(i)%row_indices
|
||||||
if (isloggedlhs(eq)) then
|
x_all => x
|
||||||
x(var) = exp(log(x(var)) - fvec(eq))
|
x_block = x(x_indices)
|
||||||
else
|
call trust_region_solve(x_block, matlab_fcn, info, tolx, tolf, maxiter, factor)
|
||||||
x(var) = x(var) - fvec(eq)
|
x(x_indices) = x_block
|
||||||
end if
|
end block
|
||||||
else
|
end do
|
||||||
if (debug) then
|
|
||||||
write (debug_msg, "('LHS variable is not determined by RHS expression (', i0, ')')") eq
|
! Verify that we have a solution
|
||||||
call mexPrintf_trim_newline(debug_msg)
|
! Note that here we use the ∞-norm, while trust region uses 2-norm; otherwise
|
||||||
write (debug_msg, "(a, ' -> ', a)") lhs_eq, endo_name_var
|
! this check would almost always fail (because the 2-norm of the full fvec is
|
||||||
call mexPrintf_trim_newline(debug_msg)
|
! larger than the 2-norm of its sub-vectors)
|
||||||
end if
|
! If the check fails, this normally means that the block decomposition was
|
||||||
if (lhs_eq(1:9) == "AUX_DIFF_" .or. lhs_eq(1:13) == "log(AUX_DIFF_") then
|
! incorrect (because some element of the Jacobian was numerically zero at the
|
||||||
if (isauxdiffloggedrhs(eq)) then
|
! guess value while not being symbolically zero)
|
||||||
x(var) = exp(log(x(var)) + fvec(eq))
|
nullify(x_indices, f_indices, x_all)
|
||||||
else
|
call matlab_fcn(x, fvec)
|
||||||
x(var) = x(var) + fvec(eq)
|
if (maxval(abs(fvec)) > tolf) then
|
||||||
end if
|
if (debug) &
|
||||||
else
|
call mexPrintf_trim_newline("DYNARE_SOLVE (solve_algo=13): residuals&
|
||||||
call mexErrMsgTxt("Algorithm solve_algo=14 cannot be used with this nonlinear problem")
|
& still too large, solving for the whole model")
|
||||||
end if
|
call trust_region_solve(x, matlab_fcn, info, tolx, tolf, maxiter, factor)
|
||||||
end if
|
|
||||||
end block
|
|
||||||
end associate
|
|
||||||
cycle
|
|
||||||
else
|
else
|
||||||
if (debug) then
|
if (size(blocks) > 1) then
|
||||||
write (debug_msg, "('DYNARE_SOLVE (solve_algo=13|14): solving block ', i0, ' with trust region routine')") i
|
! Note that the value of info may be different across blocks
|
||||||
call mexPrintf_trim_newline(debug_msg)
|
info = 1
|
||||||
end if
|
end if
|
||||||
end if
|
end if
|
||||||
|
else ! No block decomposition
|
||||||
if (size(blocks(i)%col_indices) /= size(blocks(i)%row_indices)) then
|
nullify(x_indices, f_indices, x_all)
|
||||||
! 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.")
|
|
||||||
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")
|
|
||||||
call trust_region_solve(x, matlab_fcn, info, tolx, tolf, maxiter, factor)
|
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
|
end if
|
||||||
|
|
||||||
plhs(1) = mxCreateDoubleMatrix(int(size(x, 1), mwSize), 1_mwSize, mxREAL)
|
plhs(1) = mxCreateDoubleMatrix(int(size(x, 1), mwSize), 1_mwSize, mxREAL)
|
||||||
|
|
|
||||||
|
|
@ -1 +1 @@
|
||||||
Subproject commit f725c534ef1344813421c2404de2cb4f3774d113
|
Subproject commit c48248fc0d65ed70fb13dd508ed0254f113cbcb6
|
||||||
|
|
@ -498,7 +498,19 @@ MODFILES = \
|
||||||
log_transform/example1.mod \
|
log_transform/example1.mod \
|
||||||
log_transform/ramst.mod \
|
log_transform/ramst.mod \
|
||||||
log_transform/fs2000_nonstationary.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 = \
|
ECB_MODFILES = \
|
||||||
var-expectations/1/example1.mod \
|
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.m.trs: ep/rbcii.m.trs
|
||||||
ep/rbcii_MCP.o.trs: ep/rbcii.o.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
|
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)))
|
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)))
|
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)))
|
m/model-inversion: $(patsubst %.mod, %.m.trs, $(filter model-inversion/%.mod, $(MODFILES)))
|
||||||
o/model-inversion: $(patsubst %.mod, %.o.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
|
# ECB files
|
||||||
M_ECB_TRS_FILES = $(patsubst %.mod, %.m.trs, $(ECB_MODFILES))
|
M_ECB_TRS_FILES = $(patsubst %.mod, %.m.trs, $(ECB_MODFILES))
|
||||||
|
|
||||||
|
|
@ -1448,7 +1486,12 @@ EXTRA_DIST = \
|
||||||
solver-test-functions/variablydimensioned.m \
|
solver-test-functions/variablydimensioned.m \
|
||||||
solver-test-functions/watson.m \
|
solver-test-functions/watson.m \
|
||||||
solver-test-functions/wood.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
|
if ENABLE_MATLAB
|
||||||
check-local: check-matlab
|
check-local: check-matlab
|
||||||
|
|
@ -1639,6 +1682,8 @@ clean-local:
|
||||||
|
|
||||||
rm -rf tests/pac/var-12/toto
|
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 "*.tex" -type f -delete
|
||||||
find . -name "*.aux" -type f -delete
|
find . -name "*.aux" -type f -delete
|
||||||
find . -name "*.log" -type f -delete
|
find . -name "*.log" -type f -delete
|
||||||
|
|
|
||||||
|
|
@ -31,8 +31,6 @@ tolx = 1e-6;
|
||||||
maxit = 50;
|
maxit = 50;
|
||||||
factor = 10;
|
factor = 10;
|
||||||
|
|
||||||
auxstruct = struct();
|
|
||||||
|
|
||||||
% List of function handles
|
% List of function handles
|
||||||
objfun = { @rosenbrock,
|
objfun = { @rosenbrock,
|
||||||
@powell1,
|
@powell1,
|
||||||
|
|
@ -73,7 +71,7 @@ for i=1:length(objfun)
|
||||||
x = objfun{i}();
|
x = objfun{i}();
|
||||||
end
|
end
|
||||||
try
|
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 isequal(func2str(objfun{i}), 'powell2')
|
||||||
if ~errorflag
|
if ~errorflag
|
||||||
testFailed = testFailed+1;
|
testFailed = testFailed+1;
|
||||||
|
|
|
||||||
|
|
@ -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;
|
||||||
|
|
@ -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
|
||||||
|
|
||||||
|
|
||||||
|
|
@ -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
|
||||||
|
|
||||||
|
|
||||||
|
|
@ -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);
|
||||||
|
|
@ -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;
|
||||||
|
|
@ -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
|
||||||
|
|
||||||
|
|
||||||
|
|
@ -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
|
||||||
|
|
||||||
|
|
||||||
|
|
@ -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);
|
||||||
|
|
@ -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;
|
||||||
|
|
@ -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
|
||||||
|
|
@ -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
|
||||||
|
|
@ -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);
|
||||||
|
|
@ -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;
|
||||||
|
|
@ -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
|
||||||
|
|
@ -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
|
||||||
|
|
@ -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
|
||||||
|
});
|
||||||
|
|
@ -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
|
||||||
Loading…
Reference in New Issue