From 3babb7ebfcc59879ff7f3e9207c34c455755de6d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Wed, 12 Mar 2025 16:47:59 +0100
Subject: [PATCH] WIP: rewrite with intermediate data in oo_
[skip ci]
---
.../controlled_paths_by_period.m | 44 +++++++++++++++++--
...rolled_paths_substitute_stacked_jacobian.m | 14 +++---
.../perfect_foresight_setup.m | 6 +++
.../perfect_foresight_solver.m | 43 +-----------------
.../perfect_foresight_solver_core.m | 7 +--
matlab/perfect-foresight-models/sim1.m | 37 ++++++++--------
6 files changed, 76 insertions(+), 75 deletions(-)
diff --git a/matlab/perfect-foresight-models/controlled_paths_by_period.m b/matlab/perfect-foresight-models/controlled_paths_by_period.m
index a2eb7b99f..6d72bf2fc 100644
--- a/matlab/perfect-foresight-models/controlled_paths_by_period.m
+++ b/matlab/perfect-foresight-models/controlled_paths_by_period.m
@@ -1,4 +1,4 @@
-function flips_by_period = controlled_paths_by_period(M_, options_)
+function controlled_paths_by_period = controlled_paths_by_period(M_, options_)
% Reorganize M_.perfect_foresight_controlled_paths period by period
% Copyright © 2025 Dynare Team
@@ -20,14 +20,50 @@ function flips_by_period = controlled_paths_by_period(M_, options_)
[periods, first_simulation_period] = get_simulation_periods(options_);
-flips_by_period = struct('exogenize_id', cell(periods, 1), 'endogenize_id', cell(periods, 1));
+controlled_paths_by_period = struct('exogenize_id', cell(periods, 1), 'endogenize_id', cell(periods, 1), 'values', cell(periods, 1));
for i=1:length(M_.perfect_foresight_controlled_paths)
prange = M_.perfect_foresight_controlled_paths(i).periods;
if isa(prange, 'dates')
prange = transpose(prange - first_simulation_period + 1);
end
for p = prange
- flips_by_period(p).exogenize_id(end+1) = M_.perfect_foresight_controlled_paths(i).exogenize_id;
- flips_by_period(p).endogenize_id(end+1) = M_.perfect_foresight_controlled_paths(i).endogenize_id;
+ controlled_paths_by_period(p).exogenize_id(end+1) = M_.perfect_foresight_controlled_paths(i).exogenize_id;
+ controlled_paths_by_period(p).endogenize_id(end+1) = M_.perfect_foresight_controlled_paths(i).endogenize_id;
+ controlled_paths_by_period(p).values(end+1) = M_.perfect_foresight_controlled_paths(i).value;
+ end
+end
+
+% Do various sanity checks
+for p = 1:periods
+ exogenize_id = controlled_paths_by_period(p).exogenize_id;
+ [~, idx] = unique(exogenize_id, 'stable');
+ duplicate_idx = setdiff(1:numel(exogenize_id), idx);
+ if ~isempty(duplicate_idx)
+ error('perfect_foresight_controlled_paths: variable %s is exogenized two times in period %d', M_.endo_names{exogenize_id(duplicate_idx(1))}, p);
+ end
+
+ endogenize_id = controlled_paths_by_period(p).endogenize_id;
+ [~, idx] = unique(endogenize_id, 'stable');
+ duplicate_idx = setdiff(1:numel(endogenize_id), idx);
+ if ~isempty(duplicate_idx)
+ error('perfect_foresight_controlled_paths: variable %s is endogenized two times in period %d', M_.exo_names{endogenize_id(duplicate_idx(1))}, p);
+ end
+end
+
+if isfield(M_, 'det_shocks')
+ for i = 1:length(M_.det_shocks)
+ if M_.det_shocks(i).exo_det
+ continue
+ end
+ prange = M_.det_shocks(i).periods;
+ if isa(prange, 'dates')
+ prange = transpose(prange - first_simulation_period + 1);
+ end
+ for p = prange
+ exo_id = M_.det_shocks(i).exo_id;
+ if ismember(exo_id, controlled_paths_by_period(p).endogenize_id)
+ error('Exogenous variables %s is present in both the shocks and the perfect_foresight_controlled_paths in period %d', M_.exo_names{exo_id}, p);
+ end
+ end
end
end
diff --git a/matlab/perfect-foresight-models/controlled_paths_substitute_stacked_jacobian.m b/matlab/perfect-foresight-models/controlled_paths_substitute_stacked_jacobian.m
index 04a361d1d..1be87d245 100644
--- a/matlab/perfect-foresight-models/controlled_paths_substitute_stacked_jacobian.m
+++ b/matlab/perfect-foresight-models/controlled_paths_substitute_stacked_jacobian.m
@@ -1,4 +1,4 @@
-function A = controlled_paths_substitute_stacked_jacobian(A, y, y0, yT, exogenousvariables, steadystate, M_, options_)
+function A = controlled_paths_substitute_stacked_jacobian(A, y, y0, yT, exogenousvariables, steadystate, controlled_paths_by_period, M_)
% Given a stacked Jacobian for the perfect foresight problem, modifies some columns
% according to the perfect_foresight_controlled_paths block.
@@ -19,14 +19,12 @@ function A = controlled_paths_substitute_stacked_jacobian(A, y, y0, yT, exogenou
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-periods = get_simulation_periods(options_);
-
-flips_by_period = controlled_paths_by_period(M_, options_);
-
dynamic_g1 = str2func([M_.fname '.sparse.dynamic_g1']);
+periods = numel(controlled_paths_by_period);
+
for p = 1:periods
- nsubst = length(flips_by_period(p).exogenize_id);
+ nsubst = length(controlled_paths_by_period(p).exogenize_id);
if nsubst > 0
if p == 1
y3n = [ y0; y(1:2*M_.endo_nbr) ];
@@ -38,8 +36,8 @@ for p = 1:periods
g1 = dynamic_g1(y3n, exogenousvariables(p+M_.maximum_lag,:), M_.params, steadystate, M_.dynamic_g1_sparse_rowval, M_.dynamic_g1_sparse_colval, M_.dynamic_g1_sparse_colptr);
subst = [ zeros((p-1)*M_.endo_nbr, nsubst);
- g1(:,3*M_.endo_nbr+flips_by_period(p).endogenize_id);
+ g1(:,3*M_.endo_nbr+controlled_paths_by_period(p).endogenize_id);
zeros((periods-p)*M_.endo_nbr, nsubst) ];
- A(:,(p-1)*M_.endo_nbr+flips_by_period(p).exogenize_id) = subst;
+ A(:,(p-1)*M_.endo_nbr+controlled_paths_by_period(p).exogenize_id) = subst;
end
end
diff --git a/matlab/perfect-foresight-models/perfect_foresight_setup.m b/matlab/perfect-foresight-models/perfect_foresight_setup.m
index ff65afb79..1ac6e56e0 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_setup.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_setup.m
@@ -74,3 +74,9 @@ end
oo_ = make_ex_(M_,options_,oo_);
oo_ = make_y_(M_,options_,oo_);
+
+if isempty(M_.perfect_foresight_controlled_paths)
+ oo_.controlled_paths_by_period = [];
+else
+ oo_.controlled_paths_by_period = controlled_paths_by_period(M_, options_);
+end
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver.m b/matlab/perfect-foresight-models/perfect_foresight_solver.m
index 329c08605..1936c08b4 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver.m
@@ -107,45 +107,6 @@ if options_.simul.homotopy_max_completion_share < 1 && ~options_.simul.homotopy_
error('perfect_foresight_solver: Option homotopy_max_completion_share has a value less than 1, so you must also specify either homotopy_linearization_fallback or homotopy_marginal_linearization_fallback')
end
-if ~isempty(M_.perfect_foresight_controlled_paths)
- flips_by_period = controlled_paths_by_period(M_, options_);
-
- for p = 1:periods
- exogenize_id = flips_by_period(p).exogenize_id;
- [~, idx] = unique(exogenize_id, 'stable');
- duplicate_idx = setdiff(1:numel(exogenize_id), idx);
- if ~isempty(duplicate_idx)
- error('perfect_foresight_controlled_paths: variable %s is exogenized two times in period %d', M_.endo_names{exogenize_id(duplicate_idx(1))}, p);
- end
-
- endogenize_id = flips_by_period(p).endogenize_id;
- [~, idx] = unique(endogenize_id, 'stable');
- duplicate_idx = setdiff(1:numel(endogenize_id), idx);
- if ~isempty(duplicate_idx)
- error('perfect_foresight_controlled_paths: variable %s is endogenized two times in period %d', M_.exo_names{endogenize_id(duplicate_idx(1))}, p);
- end
- end
-
- if isfield(M_, 'det_shocks')
- for i = 1:length(M_.det_shocks)
- if M_.det_shocks(i).exo_det
- continue
- end
- prange = M_.det_shocks(i).periods;
- if isa(prange, 'dates')
- prange = transpose(prange - first_simulation_period + 1);
- end
- for p = prange
- exo_id = M_.det_shocks(i).exo_id;
- if ismember(exo_id, flips_by_period(p).endogenize_id)
- error('Exogenous variables %s is present in both the shocks and the perfect_foresight_controlled_paths in period %d', M_.exo_names{exo_id}, p);
- end
- end
- end
- end
-end
-
-
initperiods = 1:M_.maximum_lag;
simperiods = M_.maximum_lag+(1:periods);
lastperiods = M_.maximum_lag+periods+(1:M_.maximum_lead);
@@ -260,7 +221,7 @@ elseif options_.simul.homotopy_marginal_linearization_fallback > 0 && completed_
extra_simul_time_counter = tic;
[extra_success, extra_endo_simul, extra_exo_simul, extra_steady_state, extra_exo_steady_state] = create_scenario(M_,options_,oo_,extra_share, shareorig, endoorig, exoorig, endobase, exobase, initperiods, lastperiods, recompute_final_steady_state, endo_simul, steady_state, exo_steady_state);
if extra_success
- [extra_endo_simul, extra_success] = perfect_foresight_solver_core(extra_endo_simul, extra_exo_simul, extra_steady_state, extra_exo_steady_state, M_, options_);
+ [extra_endo_simul, extra_success] = perfect_foresight_solver_core(extra_endo_simul, extra_exo_simul, extra_steady_state, extra_exo_steady_state, oo_.controlled_paths_by_period, M_, options_);
end
if ~extra_success
if ~options_.noprint
@@ -392,7 +353,7 @@ while step > options_.simul.homotopy_min_step_size
end
% Solve for the paths of the endogenous variables.
- [endo_simul, success, maxerror, solver_iter, per_block_status, exo_simul] = perfect_foresight_solver_core(endo_simul, exo_simul, steady_state, exo_steady_state, M_, options_);
+ [endo_simul, success, maxerror, solver_iter, per_block_status, exo_simul] = perfect_foresight_solver_core(endo_simul, exo_simul, steady_state, exo_steady_state, oo_.controlled_paths_by_period, M_, options_);
else
success = false;
maxerror = NaN;
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
index 8d38504c3..82795a05b 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
@@ -1,4 +1,4 @@
-function [y, success, maxerror, iter, per_block_status, exo_simul] = perfect_foresight_solver_core(y, exo_simul, steady_state, exo_steady_state, M_, options_)
+function [y, success, maxerror, iter, per_block_status, exo_simul] = perfect_foresight_solver_core(y, exo_simul, steady_state, exo_steady_state, controlled_paths_by_period, M_, options_)
% Core function calling solvers for perfect foresight model
%
@@ -7,6 +7,7 @@ function [y, success, maxerror, iter, per_block_status, exo_simul] = perfect_for
% - exo_simul [matrix] path of exogenous
% - steady_state [vector] steady state of endogenous variables
% - exo_steady_state [vector] steady state of exogenous variables
+% - controlled_paths_by_period [struct] data from perfect_foresight_controlled_paths block
% - M_ [struct] contains a description of the model.
% - options_ [struct] contains various options.
%
@@ -55,7 +56,7 @@ if options_.linear && (isequal(options_.stack_solve_algo, 0) || isequal(options_
options_.linear_approximation = true;
end
-if ~isempty(M_.perfect_foresight_controlled_paths)
+if ~isempty(controlled_paths_by_period)
assert(nargout >= 6); % Ensure modified exos are used
if options_.bytecode
error('perfect_foresight_controlled_paths is not compatible with bytecode option')
@@ -122,7 +123,7 @@ else
end
[y, success, maxerror] = sim1_linear(y, exo_simul, steady_state, exo_steady_state, M_, options_);
else
- [y, success, maxerror, iter, exo_simul] = sim1(y, exo_simul, steady_state, M_, options_);
+ [y, success, maxerror, iter, exo_simul] = sim1(y, exo_simul, steady_state, controlled_paths_by_period, M_, options_);
end
case {1 6}
if options_.linear_approximation
diff --git a/matlab/perfect-foresight-models/sim1.m b/matlab/perfect-foresight-models/sim1.m
index 8b013b061..274615b02 100644
--- a/matlab/perfect-foresight-models/sim1.m
+++ b/matlab/perfect-foresight-models/sim1.m
@@ -1,4 +1,4 @@
-function [endogenousvariables, success, err, iter, exogenousvariables] = sim1(endogenousvariables, exogenousvariables, steadystate, M_, options_)
+function [endogenousvariables, success, err, iter, exogenousvariables] = sim1(endogenousvariables, exogenousvariables, steadystate, controlled_paths_by_period, M_, options_)
% [endogenousvariables, success, err, iter] = sim1(endogenousvariables, exogenousvariables, steadystate, M_, options_)
% Performs deterministic simulations with lead or lag of one period, using
% a basic Newton solver on sparse matrices.
@@ -8,6 +8,7 @@ function [endogenousvariables, success, err, iter, exogenousvariables] = sim1(en
% - 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_.maximum_lag+M_.maximum_lead)*M array, paths for the exogenous variables.
% - steadystate [double] N*1 array, steady state for the endogenous variables.
+% - controlled_paths_by_period [struct] data from perfect_foresight_controlled_paths block
% - M_ [struct] contains a description of the model.
% - options_ [struct] contains various options.
% OUTPUTS
@@ -38,16 +39,15 @@ function [endogenousvariables, success, err, iter, exogenousvariables] = sim1(en
verbose = options_.verbosity && ~options_.noprint;
ny = M_.endo_nbr;
-[periods, first_simulation_period] = get_simulation_periods(options_);
+periods = get_simulation_periods(options_);
vperiods = periods*ones(1,options_.simul.maxit);
-if ~isempty(M_.perfect_foresight_controlled_paths)
- for i=1:length(M_.perfect_foresight_controlled_paths)
- prange = M_.perfect_foresight_controlled_paths(i).periods;
- if isa(prange, 'dates')
- prange = prange - first_simulation_period + 1;
+if ~isempty(controlled_paths_by_period)
+ for p = 1:periods
+ if isempty(controlled_paths_by_period(p).exogenize_id)
+ continue
end
- endogenousvariables(M_.perfect_foresight_controlled_paths(i).exogenize_id, prange+M_.maximum_lag) = M_.perfect_foresight_controlled_paths(i).value;
+ endogenousvariables(controlled_paths_by_period(p).exogenize_id, p+M_.maximum_lag) = controlled_paths_by_period(p).values;
end
if options_.debug
@@ -90,8 +90,8 @@ for iter = 1:options_.simul.maxit
% A is the stacked Jacobian with period x equations alongs the rows and
% periods times variables (in declaration order) along the columns
- if ~isempty(M_.perfect_foresight_controlled_paths)
- A = controlled_paths_substitute_stacked_jacobian(A, y, y0, yT, exogenousvariables, steadystate, M_, options_);
+ if ~isempty(controlled_paths_by_period)
+ A = controlled_paths_substitute_stacked_jacobian(A, y, y0, yT, exogenousvariables, steadystate, controlled_paths_by_period, M_);
end
if options_.debug && iter==1
@@ -161,16 +161,15 @@ for iter = 1:options_.simul.maxit
end
y = y + dy;
- if ~isempty(M_.perfect_foresight_controlled_paths)
- for i = 1:length(M_.perfect_foresight_controlled_paths)
- prange = M_.perfect_foresight_controlled_paths(i).periods;
- if isa(prange, 'dates')
- prange = prange - first_simulation_period + 1;
+ if ~isempty(controlled_paths_by_period)
+ for p = 1:periods
+ endogenize_id = controlled_paths_by_period(p).endogenize_id;
+ exogenize_id = controlled_paths_by_period(p).exogenize_id;
+ if isempty(endogenize_id)
+ continue
end
- endogenize_id = M_.perfect_foresight_controlled_paths(i).endogenize_id;
- exogenize_id = M_.perfect_foresight_controlled_paths(i).exogenize_id;
- y(exogenize_id+(prange-1)*M_.endo_nbr) = M_.perfect_foresight_controlled_paths(i).value;
- exogenousvariables(prange+M_.maximum_lag,endogenize_id) = exogenousvariables(prange+M_.maximum_lag,endogenize_id) + dy(exogenize_id+(prange-1)*M_.endo_nbr);
+ y(exogenize_id+(p-1)*M_.endo_nbr) = controlled_paths_by_period(p).values;
+ exogenousvariables(p+M_.maximum_lag,endogenize_id) = exogenousvariables(p+M_.maximum_lag,endogenize_id) + dy(exogenize_id+(p-1)*M_.endo_nbr)';
end
end
end
--
GitLab