Verified Commit 5b591fac authored by Sébastien Villemot's avatar Sébastien Villemot
Browse files

New perfect_foresight_problem MEX file

It constructs the stacked residuals and jacobian of the perfect foresight
problem.

It is an almost perfect replacement for the perfect_foresight_problem.m
routine, while being much more efficient.

Note however that the DLL never return complex numbers (it instead puts NaNs at
the place where there would have been complex). This may create problems for
some MOD files; the algorithms will need to be adapted to use a more
line-search method.
parent 69229b6b
function [residuals,JJacobian] = perfect_foresight_problem(y, dynamic_function, Y0, YT, ...
exo_simul, params, steady_state, ...
maximum_lag, T, ny, i_cols, ...
i_cols_J1, i_cols_1, i_cols_T, ...
i_cols_j, i_cols_0, i_cols_J0)
% Computes the residuals and the Jacobian matrix for a perfect foresight problem over T periods.
%
% INPUTS
% - y [double] N*1 array, terminal conditions for the endogenous variables
% - dynamic_function [handle] function handle to _dynamic-file
% - Y0 [double] N*1 array, initial conditions for the endogenous variables
% - YT [double] N*1 array, terminal conditions for the endogenous variables
% - exo_simul [double] nperiods*M_.exo_nbr matrix of exogenous variables (in declaration order)
% for all simulation periods
% - params [double] nparams*1 array, parameter values
% - steady_state [double] endo_nbr*1 vector of steady state values
% - maximum_lag [scalar] maximum lag present in the model
% - T [scalar] number of simulation periods
% - ny [scalar] number of endogenous variables
% - i_cols [double] indices of variables appearing in M.lead_lag_incidence
% and that need to be passed to _dynamic-file
% - i_cols_J1 [double] indices of contemporaneous and forward looking variables
% appearing in M.lead_lag_incidence
% - i_cols_1 [double] indices of contemporaneous and forward looking variables in
% M.lead_lag_incidence in dynamic Jacobian (relevant in first period)
% - i_cols_T [double] columns of dynamic Jacobian related to contemporaneous and backward-looking
% variables (relevant in last period)
% - i_cols_j [double] indices of contemporaneous variables in M.lead_lag_incidence
% in dynamic Jacobian (relevant in intermediate periods)
% - i_cols_0 [double] indices of contemporaneous variables in M.lead_lag_incidence in dynamic
% Jacobian (relevant in problems with periods=1)
% - i_cols_J0 [double] indices of contemporaneous variables appearing in M.lead_lag_incidence (relevant in problems with periods=1)
%
% OUTPUTS
% - residuals [double] (N*T)*1 array, residuals of the stacked problem
% - JJacobian [double] (N*T)*(N*T) array, Jacobian of the stacked problem
%
% ALGORITHM
% None
%
% SPECIAL REQUIREMENTS
% None.
% Copyright (C) 1996-2019 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 <http://www.gnu.org/licenses/>.
YY = [Y0; y; YT];
residuals = zeros(T*ny,1);
if nargout == 2
iJacobian = cell(T,1);
end
i_rows = 1:ny;
i_cols_J = i_cols;
offset = 0;
for it = maximum_lag+(1:T)
if nargout == 1
residuals(i_rows) = dynamic_function(YY(i_cols), exo_simul, params, steady_state, it);
elseif nargout == 2
[residuals(i_rows),jacobian] = dynamic_function(YY(i_cols), exo_simul, params, steady_state, it);
if T==1 && it==maximum_lag+1
[rows, cols, vals] = find(jacobian(:,i_cols_0));
iJacobian{1} = [rows, i_cols_J0(cols), vals];
elseif it == maximum_lag+1
[rows,cols,vals] = find(jacobian(:,i_cols_1));
iJacobian{1} = [offset+rows, i_cols_J1(cols), vals];
elseif it == maximum_lag+T
[rows,cols,vals] = find(jacobian(:,i_cols_T));
iJacobian{T} = [offset+rows, i_cols_J(i_cols_T(cols)), vals];
else
[rows,cols,vals] = find(jacobian(:,i_cols_j));
iJacobian{it-maximum_lag} = [offset+rows, i_cols_J(cols), vals];
i_cols_J = i_cols_J + ny;
end
offset = offset + ny;
end
i_rows = i_rows + ny;
i_cols = i_cols + ny;
end
if nargout == 2
iJacobian = cat(1,iJacobian{:});
JJacobian = sparse(iJacobian(:,1), iJacobian(:,2), iJacobian(:,3), T*ny, T*ny);
end
......@@ -176,26 +176,23 @@ if ~oo_.deterministic_simulation.status && ~options_.no_homotopy
end
if ~isreal(oo_.endo_simul(:)) % can only happen without bytecode
y0 = real(oo_.endo_simul(:,1));
yT = real(oo_.endo_simul(:,periods+2));
yy = real(oo_.endo_simul(:,2:periods+1));
illi = M_.lead_lag_incidence';
[i_cols,~,i_cols_j] = find(illi(:));
illi = illi(:,2:3);
[i_cols_J1,~,i_cols_1] = find(illi(:));
i_cols_T = nonzeros(M_.lead_lag_incidence(1:2,:)');
if periods==1
i_cols_0 = nonzeros(M_.lead_lag_incidence(2,:)');
i_cols_J0 = find(M_.lead_lag_incidence(2,:)');
if ~isreal(oo_.endo_simul(:)) % cannot happen with bytecode or the perfect_foresight_problem DLL
if M_.maximum_lag > 0
y0 = real(oo_.endo_simul(:, M_.maximum_lag));
else
i_cols_0 = [];
i_cols_J0 = [];
y0 = NaN(ny, 1);
end
residuals = perfect_foresight_problem(yy(:),str2func([M_.fname '.dynamic']), y0, yT, ...
oo_.exo_simul,M_.params,oo_.steady_state, ...
M_.maximum_lag, periods, M_.endo_nbr, i_cols, ...
i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0);
if M_.maximum_lead > 0
yT = real(oo_.endo_simul(:, M_.maximum_lag+periods+1));
else
yT = NaN(ny, 1);
end
yy = real(oo_.endo_simul(:,M_.maximum_lag+(1:periods)));
model_dynamic_g1_nz = str2func([M_.fname,'.dynamic_g1_nz']);
[nzij_pred, nzij_current, nzij_fwrd] = model_dynamic_g1_nz();
residuals = perfect_foresight_problem(yy(:), M_.fname, sum(M_.dynamic_tmp_nbr(1:2)), y0, yT, oo_.exo_simul, M_.params, oo_.steady_state, periods, M_.endo_nbr, M_.maximum_lag, M_.maximum_endo_lag, M_.lead_lag_incidence, nzij_pred, nzij_current, nzij_fwrd, M_.has_external_function, options_.use_dll);
if max(abs(residuals))< options_.dynatol.f
oo_.deterministic_simulation.status = 1;
oo_.endo_simul=real(oo_.endo_simul);
......
......@@ -120,31 +120,27 @@ else
end
if nargout>1
y0 = oo_.endo_simul(:,1);
yT = oo_.endo_simul(:,periods+2);
yy = oo_.endo_simul(:,2:periods+1);
illi = M_.lead_lag_incidence';
[i_cols,~,i_cols_j] = find(illi(:));
illi = illi(:,2:3);
[i_cols_J1,~,i_cols_1] = find(illi(:));
i_cols_T = nonzeros(M_.lead_lag_incidence(1:2,:)');
if periods==1
i_cols_0 = nonzeros(M_.lead_lag_incidence(2,:)');
i_cols_J0 = find(M_.lead_lag_incidence(2,:)');
else
i_cols_0 = [];
i_cols_J0 = [];
end
if options_.block && ~options_.bytecode
maxerror = oo_.deterministic_simulation.error;
else
if options_.bytecode
[chck, residuals, ~]= bytecode('dynamic','evaluate', oo_.endo_simul, oo_.exo_simul, M_.params, oo_.steady_state, 1);
else
residuals = perfect_foresight_problem(yy(:),str2func([M_.fname '.dynamic']), y0, yT, ...
oo_.exo_simul,M_.params,oo_.steady_state, ...
M_.maximum_lag, periods,M_.endo_nbr,i_cols, ...
i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0);
if M_.maximum_lag > 0
y0 = oo_.endo_simul(:, M_.maximum_lag);
else
y0 = NaN(ny, 1);
end
if M_.maximum_lead > 0
yT = oo_.endo_simul(:, M_.maximum_lag+periods+1);
else
yT = NaN(ny, 1);
end
yy = oo_.endo_simul(:,M_.maximum_lag+(1:periods));
model_dynamic_g1_nz = str2func([M_.fname,'.dynamic_g1_nz']);
[nzij_pred, nzij_current, nzij_fwrd] = model_dynamic_g1_nz();
residuals = perfect_foresight_problem(yy(:), M_.fname, sum(M_.dynamic_tmp_nbr(1:2)), y0, yT, oo_.exo_simul, M_.params, oo_.steady_state, periods, M_.endo_nbr, M_.maximum_lag, M_.maximum_endo_lag, M_.lead_lag_incidence, nzij_pred, nzij_current, nzij_fwrd, M_.has_external_function, options_.use_dll);
end
maxerror = max(max(abs(residuals)));
end
......
......@@ -70,8 +70,16 @@ elseif (options.solve_algo == 11)
end
end
y0 = endogenousvariables(:,M.maximum_lag);
yT = endogenousvariables(:,M.maximum_lag+periods+1);
if M.maximum_lag > 0
y0 = endogenousvariables(:, M.maximum_lag);
else
y0 = NaN(ny, 1);
end
if M.maximum_lead > 0
yT = endogenousvariables(:, M.maximum_lag+periods+1);
else
yT = NaN(ny, 1);
end
z = endogenousvariables(:,M.maximum_lag+(1:periods));
illi = M.lead_lag_incidence';
[i_cols,~,i_cols_j] = find(illi(:));
......@@ -85,4 +93,4 @@ else
i_cols_0 = [];
i_cols_J0 = [];
end
dynamicmodel = str2func([M.fname,'.dynamic']);
\ No newline at end of file
dynamicmodel = str2func([M.fname,'.dynamic']);
function [endogenousvariables, info] = sim1(endogenousvariables, exogenousvariables, steadystate, M, options)
% Performs deterministic simulations with lead or lag on one period. Uses sparse matrices.
% Performs deterministic simulations with lead or lag of one period, using
% a basic Newton solver on sparse matrices.
% Uses perfect_foresight_problem DLL to construct the stacked problem.
%
% INPUTS
% - endogenousvariables [double] N*T array, paths for the endogenous variables (initial guess).
% - endogenousvariables [double] N*(T+M.maximum_lag+M.maximum_lead) array, paths for the endogenous variables (initial condition + initial guess + terminal condition).
% - exogenousvariables [double] T*M array, paths for the exogenous variables.
% - steadystate [double] N*1 array, steady state for the endogenous variables.
% - M [struct] contains a description of the model.
% - options [struct] contains various options.
% OUTPUTS
% - endogenousvariables [double] N*T array, paths for the endogenous variables (solution of the perfect foresight model).
% - endogenousvariables [double] N*(T+M.maximum_lag+M.maximum_lead) array, paths for the endogenous variables (solution of the perfect foresight model).
% - info [struct] contains informations about the results.
% ALGORITHM
% ...
%
% SPECIAL REQUIREMENTS
% None.
% Copyright (C) 1996-2019 Dynare Team
%
......@@ -36,38 +32,25 @@ function [endogenousvariables, info] = sim1(endogenousvariables, exogenousvariab
verbose = options.verbosity;
endogenous_terminal_period = options.endogenous_terminal_period;
vperiods = options.periods*ones(1,options.simul.maxit);
azero = options.dynatol.f/1e7;
lead_lag_incidence = M.lead_lag_incidence;
ny = M.endo_nbr;
periods = options.periods;
vperiods = periods*ones(1,options.simul.maxit);
maximum_lag = M.maximum_lag;
max_lag = M.maximum_endo_lag;
nyp = nnz(lead_lag_incidence(1,:)) ;
ny0 = nnz(lead_lag_incidence(2,:)) ;
nyf = nnz(lead_lag_incidence(3,:)) ;
nd = nyp+ny0+nyf;
stop = 0 ;
if M.maximum_lag > 0
y0 = endogenousvariables(:, M.maximum_lag);
else
y0 = NaN(ny, 1);
end
periods = options.periods;
params = M.params;
if M.maximum_lead > 0
yT = endogenousvariables(:, M.maximum_lag+periods+1);
else
yT = NaN(ny, 1);
end
i_cols_1 = nonzeros(lead_lag_incidence(2:3,:)');
i_cols_A1 = find(lead_lag_incidence(2:3,:)');
i_cols_A1 = i_cols_A1(:);
i_cols_T = nonzeros(lead_lag_incidence(1:2,:)');
i_cols_0 = nonzeros(lead_lag_incidence(2,:)');
i_cols_A0 = find(lead_lag_incidence(2,:)');
i_cols_A0 = i_cols_A0(:);
i_cols_j = (1:nd)';
i_upd = maximum_lag*ny+(1:periods*ny);
y = reshape(endogenousvariables(:, M.maximum_lag+(1:periods)), ny*periods, 1);
Y = endogenousvariables(:);
stop = false;
if verbose
skipline()
......@@ -76,121 +59,78 @@ if verbose
skipline()
end
model_dynamic = str2func([M.fname,'.dynamic']);
res = zeros(periods*ny,1);
o_periods = periods;
if endogenous_terminal_period
ZERO = zeros(length(i_upd),1);
end
model_dynamic_g1_nz = str2func([M.fname,'.dynamic_g1_nz']);
[nzij_pred, nzij_current, nzij_fwrd] = model_dynamic_g1_nz();
h1 = clock ;
iA = zeros(periods*M.NNZDerivatives(1),3);
h1 = clock;
for iter = 1:options.simul.maxit
h2 = clock ;
i_rows = (1:ny)';
i_cols_A = find(lead_lag_incidence');
i_cols_A = i_cols_A(:);
i_cols = i_cols_A+(maximum_lag-1)*ny;
m = 0;
for it = (maximum_lag+1):(maximum_lag+periods)
[d1,jacobian] = model_dynamic(Y(i_cols), exogenousvariables, params, steadystate,it);
if periods==1 && it==maximum_lag+1
[r,c,v] = find(jacobian(:,i_cols_0));
iA((1:length(v))+m,:) = [i_rows(r(:)),i_cols_A0(c(:)),v(:)];
elseif it == maximum_lag+periods
[r,c,v] = find(jacobian(:,i_cols_T));
iA((1:length(v))+m,:) = [i_rows(r(:)),i_cols_A(i_cols_T(c(:))),v(:)];
elseif it == maximum_lag+1
[r,c,v] = find(jacobian(:,i_cols_1));
iA((1:length(v))+m,:) = [i_rows(r(:)),i_cols_A1(c(:)),v(:)];
else
[r,c,v] = find(jacobian(:,i_cols_j));
iA((1:length(v))+m,:) = [i_rows(r(:)),i_cols_A(c(:)),v(:)];
end
m = m + length(v);
res(i_rows) = d1;
if endogenous_terminal_period && iter>1
dr = max(abs(d1));
if dr<azero
vperiods(iter) = it;
periods = it-maximum_lag+1;
h2 = clock;
[res, A] = perfect_foresight_problem(y, M.fname, sum(M.dynamic_tmp_nbr(1:2)), y0, yT, exogenousvariables, M.params, steadystate, periods, ny, M.maximum_lag, M.maximum_endo_lag, M.lead_lag_incidence, nzij_pred, nzij_current, nzij_fwrd, M.has_external_function, options.use_dll);
if options.endogenous_terminal_period && iter > 1
for it = 1:periods
if max(abs(res((it-1)*ny+(1:ny)))) < options.dynatol.f/1e7
if it < periods
res = res(1:(it*ny));
A = A(1:(it*ny), 1:(it*ny));
yT = y(it*ny+(1:ny));
endogenousvariables(:, M.maximum_lag+((it+1):periods)) = reshape(y(it*ny+(1:(ny*(periods-it)))), ny, periods-it);
y = y(1:(it*ny));
periods = it;
end
break
end
end
i_rows = i_rows + ny;
i_cols = i_cols + ny;
if it > maximum_lag+1
i_cols_A = i_cols_A + ny;
end
vperiods(iter) = periods;
end
err = max(abs(res));
if options.debug
fprintf('\nLargest absolute residual at iteration %d: %10.3f\n',iter,err);
if any(isnan(res)) || any(isinf(res)) || any(isnan(Y)) || any(isinf(Y))
if any(isnan(res)) || any(isinf(res)) || any(any(isnan(endogenousvariables))) || any(any(isinf(endogenousvariables)))
fprintf('\nWARNING: NaN or Inf detected in the residuals or endogenous variables.\n');
end
if ~isreal(res) || ~isreal(Y)
fprintf('\nWARNING: Imaginary parts detected in the residuals or endogenous variables.\n');
end
skipline()
end
if verbose
str = sprintf('Iter: %s,\t err. = %s, \t time = %s',num2str(iter),num2str(err), num2str(etime(clock,h2)));
disp(str);
fprintf('Iter: %d,\t err. = %g,\t time = %g\n', iter, err, etime(clock,h2));
end
if err < options.dynatol.f
stop = 1 ;
stop = true;
break
end
iA = iA(1:m,:);
A = sparse(iA(:,1),iA(:,2),iA(:,3),periods*ny,periods*ny);
if endogenous_terminal_period && iter>1
dy = ZERO;
if options.simul.robust_lin_solve
dy(1:i_rows(end)) = -lin_solve_robust( A(1:i_rows(end),1:i_rows(end)), res(1:i_rows(end)),verbose );
else
dy(1:i_rows(end)) = -lin_solve( A(1:i_rows(end),1:i_rows(end)), res(1:i_rows(end)), verbose );
end
if options.simul.robust_lin_solve
dy = -lin_solve_robust(A, res, verbose);
else
if options.simul.robust_lin_solve
dy = -lin_solve_robust( A, res, verbose );
else
dy = -lin_solve( A, res, verbose );
end
dy = -lin_solve(A, res, verbose);
end
if any(~isreal(dy)) || any(isnan(dy)) || any(isinf(dy))
if any(isnan(dy)) || any(isinf(dy))
if verbose
display_critical_variables(reshape(dy,[ny periods])', M);
end
end
Y(i_upd) = Y(i_upd) + dy;
y = y + dy;
end
if endogenous_terminal_period
err = evaluate_max_dynamic_residual(model_dynamic, Y, exogenousvariables, params, steadystate, o_periods, ny, max_lag, lead_lag_incidence);
periods = o_periods;
end
endogenousvariables(:, M.maximum_lag+(1:periods)) = reshape(y, ny, periods);
if options.endogenous_terminal_period
periods = options.periods;
err = evaluate_max_dynamic_residual(str2func([M.fname,'.dynamic']), endogenousvariables, exogenousvariables, M.params, steadystate, periods, ny, M.maximum_endo_lag, M.lead_lag_incidence);
end
if stop
if any(isnan(res)) || any(isinf(res)) || any(isnan(Y)) || any(isinf(Y)) || ~isreal(res) || ~isreal(Y)
if any(any(isnan(endogenousvariables))) || any(any(isinf(endogenousvariables)))
info.status = false;% NaN or Inf occurred
info.error = err;
info.iterations = iter;
info.periods = vperiods(1:iter);
endogenousvariables = reshape(Y,ny,periods+maximum_lag+M.maximum_lead);
if verbose
skipline()
disp(sprintf('Total time of simulation: %s.', num2str(etime(clock,h1))))
if ~isreal(res) || ~isreal(Y)
disp('Simulation terminated with imaginary parts in the residuals or endogenous variables.')
else
disp('Simulation terminated with NaN or Inf in the residuals or endogenous variables.')
end
fprintf('Total time of simulation: %g.\n', etime(clock,h1))
disp('Simulation terminated with NaN or Inf in the residuals or endogenous variables.')
display_critical_variables(reshape(dy,[ny periods])', M);
disp('There is most likely something wrong with your model. Try model_diagnostics or another simulation method.')
printline(105)
......@@ -198,19 +138,18 @@ if stop
else
if verbose
skipline();
disp(sprintf('Total time of simulation: %s', num2str(etime(clock,h1))))
fprintf('Total time of simulation: %g.\n', etime(clock,h1))
printline(56)
end
info.status = true;% Convergency obtained.
info.error = err;
info.iterations = iter;
info.periods = vperiods(1:iter);
endogenousvariables = reshape(Y,ny,periods+maximum_lag+M.maximum_lead);
end
elseif ~stop
if verbose
skipline();
disp(sprintf('Total time of simulation: %s.', num2str(etime(clock,h1))))
fprintf('Total time of simulation: %g.\n', etime(clock,h1))
disp('Maximum number of iterations is reached (modify option maxit).')
printline(62)
end
......@@ -224,21 +163,21 @@ if verbose
skipline();
end
function x = lin_solve( A, b,verbose)
if norm( b ) < sqrt( eps ) % then x = 0 is a solution
function x = lin_solve(A, b, verbose)
if norm(b) < sqrt(eps) % then x = 0 is a solution
x = 0;
return
end
x = A\b;
x( ~isfinite( x ) ) = 0;
relres = norm( b - A * x ) / norm( b );
x(~isfinite(x)) = 0;
relres = norm(b - A*x) / norm(b);
if relres > 1e-6 && verbose
fprintf( 'WARNING : Failed to find a solution to the linear system.\n' );
fprintf('WARNING : Failed to find a solution to the linear system.\n');
end
function [ x, flag, relres ] = lin_solve_robust( A, b , verbose)
if norm( b ) < sqrt( eps ) % then x = 0 is a solution
function [ x, flag, relres ] = lin_solve_robust(A, b ,verbose)
if norm(b) < sqrt(eps) % then x = 0 is a solution
x = 0;
flag = 0;
relres = 0;
......@@ -246,8 +185,8 @@ if norm( b ) < sqrt( eps ) % then x = 0 is a solution
end
x = A\b;
x( ~isfinite( x ) ) = 0;
[ x, flag, relres ] = bicgstab( A, b, [], [], [], [], x ); % returns immediately if x is a solution
x(~isfinite(x)) = 0;
[ x, flag, relres ] = bicgstab(A, b, [], [], [], [], x); % returns immediately if x is a solution
if flag == 0
return
end
......@@ -255,38 +194,38 @@ end
disp( relres );
if verbose
fprintf( 'Initial bicgstab failed, trying alternative start point.\n' );
fprintf('Initial bicgstab failed, trying alternative start point.\n');
end
old_x = x;
old_relres = relres;
[ x, flag, relres ] = bicgstab( A, b );
[ x, flag, relres ] = bicgstab(A, b);
if flag == 0
return
end
if verbose
fprintf( 'Alternative start point also failed with bicgstab, trying gmres.\n' );
fprintf('Alternative start point also failed with bicgstab, trying gmres.\n');
end
if old_relres < relres
x = old_x;
end
[ x, flag, relres ] = gmres( A, b, [], [], [], [], [], x );
[ x, flag, relres ] = gmres(A, b, [], [], [], [], [], x);
if flag == 0
return
end
if verbose
fprintf( 'Initial gmres failed, trying alternative start point.\n' );
fprintf('Initial gmres failed, trying alternative start point.\n');
end
old_x = x;
old_relres = relres;
[ x, flag, relres ] = gmres( A, b );
[ x, flag, relres ] = gmres(A, b);
if flag == 0
return
end
if verbose
fprintf( 'Alternative start point also failed with gmres, using the (SLOW) Moore-Penrose Pseudo-Inverse.\n' );
fprintf('Alternative start point also failed with gmres, using the (SLOW) Moore-Penrose Pseudo-Inverse.\n');
end
if old_relres < relres
x = old_x;
......@@ -294,15 +233,15 @@ if old_relres < relres
end
old_x = x;
old_relres = relres;
x = pinv( full( A ) ) * b;
relres = norm( b - A * x ) / norm( b );
x = pinv(full(A)) * b;
relres = norm(b - A*x) / norm(b);
if old_relres < relres