From 5ccfdd54003331cf680b11d68ccb6cb33d3d35fc Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Thu, 26 Sep 2024 18:16:53 +0200
Subject: [PATCH] Linear perfect foresight solver: use sparse representation
for the model
Ref. #1859
---
.../linear_perfect_foresight_problem.m | 87 ++++++-------------
.../private/initialize_stacked_problem.m | 78 -----------------
.../solve_stacked_linear_problem.m | 38 ++++----
3 files changed, 45 insertions(+), 158 deletions(-)
delete mode 100644 matlab/perfect-foresight-models/private/initialize_stacked_problem.m
diff --git a/matlab/perfect-foresight-models/linear_perfect_foresight_problem.m b/matlab/perfect-foresight-models/linear_perfect_foresight_problem.m
index adc8045aa6..439ea98c1f 100644
--- a/matlab/perfect-foresight-models/linear_perfect_foresight_problem.m
+++ b/matlab/perfect-foresight-models/linear_perfect_foresight_problem.m
@@ -1,22 +1,10 @@
function [residuals,JJacobian] = linear_perfect_foresight_problem(y, dynamicjacobian, 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, jendo, jexog)
+ exo_simul, params, steady_state, ...
+ maximum_lag, T, ny)
% Computes the residuals and the Jacobian matrix for a linear perfect foresight problem over T periods.
-%
-% INPUTS
-% ...
-%
-% OUTPUTS
-% ...
-%
-% ALGORITHM
-% ...
-%
-% SPECIAL REQUIREMENTS
-% None.
-% Copyright © 2015-2020 Dynare Team
+% Copyright © 2015-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -33,64 +21,39 @@ function [residuals,JJacobian] = linear_perfect_foresight_problem(y, dynamicjaco
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+YY = reshape([Y0; y; YT], ny, T+2);
-YY = [Y0; y; YT];
-
-residuals = zeros(T*ny,1);
-
-z = zeros(columns(dynamicjacobian), 1);
-
+residuals = NaN(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)
- z(jendo) = YY(i_cols);
- z(jexog) = transpose(exo_simul(it,:));
- residuals(i_rows) = dynamicjacobian*z;
+for t = 1:T
+ z = [ vec(YY(:, t+(0:2))); transpose(exo_simul(t+maximum_lag, :)) ];
+ residuals(offset + (1:ny)) = dynamicjacobian*z;
+
if nargout == 2
- if T==1 && it==maximum_lag+1
- [rows, cols, vals] = find(dynamicjacobian(:,i_cols_0));
- if size(dynamicjacobian, 1) == 1 % find() will return row vectors in this case
- rows = rows';
- cols = cols';
- vals = vals';
- end
- iJacobian{1} = [rows, i_cols_J0(cols), vals];
- elseif it == maximum_lag+1
- [rows,cols,vals] = find(dynamicjacobian(:,i_cols_1));
- if size(dynamicjacobian, 1) == 1 % find() will return row vectors in this case
- rows = rows';
- cols = cols';
- vals = vals';
- end
- iJacobian{1} = [offset+rows, i_cols_J1(cols), vals];
- elseif it == maximum_lag+T
- [rows,cols,vals] = find(dynamicjacobian(:,i_cols_T));
- if size(dynamicjacobian, 1) == 1 % find() will return row vectors in this case
- rows = rows';
- cols = cols';
- vals = vals';
- end
- iJacobian{T} = [offset+rows, i_cols_J(i_cols_T(cols)), vals];
+ if T == 1 % static model
+ [rows, cols, vals] = find(dynamicjacobian(:, ny+(1:ny)));
+ elseif t == 1
+ [rows, cols, vals] = find(dynamicjacobian(:, ny+(1:2*ny)));
+ elseif t == T
+ [rows, cols, vals] = find(dynamicjacobian(:, 1:2*ny));
+ cols = cols - ny;
else
- [rows,cols,vals] = find(dynamicjacobian(:,i_cols_j));
- if size(dynamicjacobian, 1) == 1 % find() will return row vectors in this case
- rows = rows';
- cols = cols';
- vals = vals';
- end
- iJacobian{it-maximum_lag} = [offset+rows, i_cols_J(cols), vals];
- i_cols_J = i_cols_J + ny;
+ [rows, cols, vals] = find(dynamicjacobian(:, 1:3*ny));
+ cols = cols - ny;
+ end
+ if size(dynamicjacobian, 1) == 1 % find() will return row vectors in this case
+ rows = rows';
+ cols = cols';
+ vals = vals';
end
- offset = offset + ny;
+ iJacobian{t} = [offset+rows, offset+cols, vals];
end
- i_rows = i_rows + ny;
- i_cols = i_cols + ny;
+ offset = offset + ny;
end
if nargout == 2
diff --git a/matlab/perfect-foresight-models/private/initialize_stacked_problem.m b/matlab/perfect-foresight-models/private/initialize_stacked_problem.m
deleted file mode 100644
index 8aea2e0f57..0000000000
--- a/matlab/perfect-foresight-models/private/initialize_stacked_problem.m
+++ /dev/null
@@ -1,78 +0,0 @@
-function [options_, y0, yT, z, i_cols, i_cols_J1, i_cols_T, i_cols_j, i_cols_1, i_cols_0, i_cols_J0, dynamicmodel] = ...
- initialize_stacked_problem(endogenousvariables, options_, M_, steadystate_y)
-
-% Sets up the stacked perfect foresight problem for use with dynare_solve.m
-%
-% INPUTS
-% - endogenousvariables [double] N*T array, paths for the endogenous variables (initial guess).
-% - options_ [struct] contains various options.
-% - M_ [struct] contains a description of the model.
-% - steadystate_y [double] N*1 array, steady state for the endogenous variables.
-%
-% OUTPUTS
-% - options_ [struct] contains various options.
-% - y0 [double] N*1 array, initial conditions for the endogenous variables
-% - yT [double] N*1 array, terminal conditions for the endogenous variables
-% - z [double] T*M array, paths for the exogenous 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_T [double] columns of dynamic Jacobian related to
-% contemporaneous and backward-looking
-% variables (relevant in last period)
-% - i_cols_j [double] indices of variables in M_.lead_lag_incidence
-% in dynamic Jacobian (relevant in intermediate periods)
-% - 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_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)
-% - dynamicmodel [handle] function handle to _dynamic-file
-
-% Copyright © 2015-2024 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/>.
-
-periods = options_.periods;
-
-if M_.maximum_lag > 0
- y0 = endogenousvariables(:, M_.maximum_lag);
-else
- y0 = NaN(M_.endo_nbr, 1);
-end
-if M_.maximum_lead > 0
- yT = endogenousvariables(:, M_.maximum_lag+periods+1);
-else
- yT = NaN(M_.endo_nbr, 1);
-end
-z = endogenousvariables(:,M_.maximum_lag+(1:periods));
-illi = M_.lead_lag_incidence';
-[i_cols,~,i_cols_j] = find(illi(:));
-if M_.maximum_endo_lag == 0
- i_cols = i_cols + M_.endo_nbr;
-end
-illi = illi(:,(1+M_.maximum_endo_lag):(1+M_.maximum_endo_lag+M_.maximum_endo_lead));
-[i_cols_J1,~,i_cols_1] = find(illi(:));
-i_cols_T = nonzeros(M_.lead_lag_incidence(1:(1+M_.maximum_endo_lag),:)');
-if periods==1
- i_cols_0 = nonzeros(M_.lead_lag_incidence(1+M_.maximum_endo_lag,:)');
- i_cols_J0 = find(M_.lead_lag_incidence(1+M_.maximum_endo_lag,:)');
-else
- i_cols_0 = [];
- i_cols_J0 = [];
-end
-dynamicmodel = str2func([M_.fname,'.dynamic']);
diff --git a/matlab/perfect-foresight-models/solve_stacked_linear_problem.m b/matlab/perfect-foresight-models/solve_stacked_linear_problem.m
index 3e104717a2..2e919a53e9 100644
--- a/matlab/perfect-foresight-models/solve_stacked_linear_problem.m
+++ b/matlab/perfect-foresight-models/solve_stacked_linear_problem.m
@@ -1,6 +1,6 @@
function [endogenousvariables, success] = solve_stacked_linear_problem(endogenousvariables, exogenousvariables, steadystate_y, steadystate_x, M_, options_)
-% Copyright © 2015-2023 Dynare Team
+% Copyright © 2015-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -17,27 +17,31 @@ function [endogenousvariables, success] = solve_stacked_linear_problem(endogenou
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-[options_, y0, yT, z, i_cols, i_cols_J1, i_cols_T, i_cols_j, i_cols_1, i_cols_0, i_cols_J0, dynamicmodel] = ...
- initialize_stacked_problem(endogenousvariables, options_, M_, steadystate_y);
-
-ip = find(M_.lead_lag_incidence(1,:)');
-ic = find(M_.lead_lag_incidence(2,:)');
-in = find(M_.lead_lag_incidence(3,:)');
+if M_.maximum_lag > 0
+ y0 = endogenousvariables(:, M_.maximum_lag);
+else
+ y0 = NaN(M_.endo_nbr, 1);
+end
+if M_.maximum_lead > 0
+ yT = endogenousvariables(:, M_.maximum_lag+options_.periods+1);
+else
+ yT = NaN(M_.endo_nbr, 1);
+end
+z = endogenousvariables(:,M_.maximum_lag+(1:options_.periods));
-% Evaluate the Jacobian of the dynamic model at the deterministic steady state.
-[d1,jacobian] = dynamicmodel(steadystate_y([ip; ic; in]), transpose(steadystate_x), M_.params, steadystate_y, 1);
+% Evaluate the residuals and Jacobian of the dynamic model at the deterministic steady state.
+y3n = repmat(steadystate_y, 3, 1);
+[d1, TT_order, TT] = feval([M_.fname,'.sparse.dynamic_resid'], y3n, steadystate_x', M_.params, ...
+ steadystate_y);
+jacobian = feval([M_.fname,'.sparse.dynamic_g1'], y3n, steadystate_x', M_.params, steadystate_y, ...
+ M_.dynamic_g1_sparse_rowval, M_.dynamic_g1_sparse_colval, ...
+ M_.dynamic_g1_sparse_colptr, TT_order, TT);
% Check that the dynamic model was evaluated at the steady state.
if ~options_.steadystate.nocheck && max(abs(d1))>1e-12
error('Jacobian is not evaluated at the steady state!')
end
-nyp = nnz(M_.lead_lag_incidence(1,:));
-ny0 = nnz(M_.lead_lag_incidence(2,:));
-nyf = nnz(M_.lead_lag_incidence(3,:));
-nd = nyp+ny0+nyf; % size of y (first argument passed to the dynamic file).
-jexog = transpose(nd+(1:M_.exo_nbr));
-jendo = transpose(1:nd);
z = bsxfun(@minus, z, steadystate_y);
x = bsxfun(@minus, exogenousvariables, steadystate_x');
@@ -46,9 +50,7 @@ x = bsxfun(@minus, exogenousvariables, steadystate_x');
options_, ...
jacobian, y0-steadystate_y, yT-steadystate_y, ...
x, M_.params, steadystate_y, ...
- M_.maximum_lag, options_.periods, M_.endo_nbr, i_cols, ...
- i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0, ...
- jendo, jexog);
+ M_.maximum_lag, options_.periods, M_.endo_nbr);
if all(imag(y)<.1*options_.dynatol.x)
if ~isreal(y)
--
GitLab