From c3d91d5ce8e93cfc02731669407776ba39c1c16d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Tue, 21 Nov 2023 18:39:33 +0100 Subject: [PATCH] Fix marginal linearization in the context of perfect_foresight_with_expectation_errors_solver with homotopy MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit If a solution corresponding to 100% of the shock can’t be found in the first informational period, marginal linearization will be performed to extrapolate a solution. However, in subsequent informational periods, this extrapolated solution cannot be used for the initial conditions of endogenous variables, because the initial conditions are not a true solution of the nonlinear model. For those subsequent informational periods, the correct approach is to compute the two solutions needed for marginal linearization using as initial conditions the values obtained in the same two solutions for the previous informational periods (stored as oo_.deterministic_simulation.{sim1,sim2}). --- .../perfect_foresight_solver.m | 61 +++++++++++++------ ...foresight_with_expectation_errors_solver.m | 13 +++- meson.build | 3 +- ...y.mod => pfwee_homotopy_linearization.mod} | 2 +- .../pfwee_homotopy_marginal_linearization.mod | 40 ++++++++++++ 5 files changed, 97 insertions(+), 22 deletions(-) rename tests/deterministic_simulations/{pfwee_homotopy.mod => pfwee_homotopy_linearization.mod} (95%) create mode 100644 tests/deterministic_simulations/pfwee_homotopy_marginal_linearization.mod diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver.m b/matlab/perfect-foresight-models/perfect_foresight_solver.m index 615dc04b6..42c8a190b 100644 --- a/matlab/perfect-foresight-models/perfect_foresight_solver.m +++ b/matlab/perfect-foresight-models/perfect_foresight_solver.m @@ -1,10 +1,15 @@ -function perfect_foresight_solver(no_error_if_learnt_in_is_present) +function perfect_foresight_solver(no_error_if_learnt_in_is_present, marginal_linearization_previous_raw_sims) % Computes deterministic simulations % % INPUTS % no_error_if_learnt_in_is_present [boolean, optional] % if true, then do not error out if a shocks(learnt_in=…) or endval(learnt_in=…) % block is present +% marginal_linearization_previous_raw_sims [struct, optional] +% if not empty, contains the two simulations used to compute the extrapolation by marginal +% linearization in a previous informational period, in the context of +% perfect_foresight_with_expectation_errors in combination with homotopy and marginal +% linearization % % OUTPUTS % none @@ -35,9 +40,12 @@ global M_ options_ oo_ check_input_arguments(options_, M_, oo_); -if nargin == 0 +if nargin < 1 no_error_if_learnt_in_is_present = false; end +if nargin < 2 + marginal_linearization_previous_raw_sims = []; +end if (~isempty(M_.learnt_shocks) || ~isempty(M_.learnt_endval)) && ~no_error_if_learnt_in_is_present error('A shocks(learnt_in=...) or endval(learnt_in=...) block is present. You want to call perfect_foresight_with_expectations_error_setup and perfect_foresight_with_expectations_error_solver.') end @@ -130,12 +138,17 @@ else recompute_final_steady_state = false; end -% Copy the paths for the exogenous and endogenous variables, as given by perfect_foresight_setup -exoorig = oo_.exo_simul; -endoorig = oo_.endo_simul; - % Perform the homotopy loop -[completed_share, endo_simul, exo_simul, steady_state, exo_steady_state, iteration, maxerror, solver_iter, per_block_status] = homotopy_loop(options_.simul.homotopy_max_completion_share, endoorig, exoorig, endobase, exobase, initperiods, simperiods, lastperiods, recompute_final_steady_state, oo_.steady_state, oo_.exo_steady_state); +if isempty(marginal_linearization_previous_raw_sims) + shareorig = 1; + endoorig = oo_.endo_simul; + exoorig = oo_.exo_simul; +else + shareorig = marginal_linearization_previous_raw_sims.sim1.homotopy_completion_share; + endoorig = marginal_linearization_previous_raw_sims.sim1.endo_simul; + exoorig = marginal_linearization_previous_raw_sims.sim1.exo_simul; +end +[completed_share, endo_simul, exo_simul, steady_state, exo_steady_state, iteration, maxerror, solver_iter, per_block_status] = homotopy_loop(options_.simul.homotopy_max_completion_share, shareorig, endoorig, exoorig, endobase, exobase, initperiods, simperiods, lastperiods, recompute_final_steady_state, oo_.steady_state, oo_.exo_steady_state); % Do linearization if needed and requested, and put results and solver status information in oo_ if completed_share == 1 @@ -180,13 +193,22 @@ elseif options_.simul.homotopy_marginal_linearization_fallback > 0 && completed_ oo_.deterministic_simulation.sim1.homotopy_completion_share = completed_share; % Now compute extra simulation. First try using the first simulation as guess value. + if isempty(marginal_linearization_previous_raw_sims) + shareorig = 1; + endoorig = oo_.endo_simul; + exoorig = oo_.exo_simul; + else + shareorig = marginal_linearization_previous_raw_sims.sim2.homotopy_completion_share; + endoorig = marginal_linearization_previous_raw_sims.sim2.endo_simul; + exoorig = marginal_linearization_previous_raw_sims.sim2.exo_simul; + end extra_share = completed_share - options_.simul.homotopy_marginal_linearization_fallback; if ~options_.noprint fprintf('Only %.1f%% of the shock could be simulated. Since marginal linearization was requested as a fallback, now running an extra simulation for %.1f%% of the shock\n\n', completed_share*100, extra_share*100) fprintf('%s\n\n', repmat('*', 1, 80)) end extra_simul_time_counter = tic; - [extra_success, extra_endo_simul, extra_exo_simul, extra_steady_state, extra_exo_steady_state] = create_scenario(extra_share, endoorig, exoorig, endobase, exobase, initperiods, lastperiods, recompute_final_steady_state, endo_simul, steady_state, exo_steady_state); + [extra_success, extra_endo_simul, extra_exo_simul, extra_steady_state, extra_exo_steady_state] = create_scenario(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_); end @@ -195,7 +217,7 @@ elseif options_.simul.homotopy_marginal_linearization_fallback > 0 && completed_ fprintf('The extra simulation for %.1f%% of the shock did not run when using the first simulation as a guess value. Now trying a full homotopy loop to get that extra simulation working\n\n', extra_share*100) fprintf('%s\n\n', repmat('*', 1, 80)) end - [extra_completed_share, extra_endo_simul, extra_exo_simul, extra_steady_state, extra_exo_steady_state] = homotopy_loop(extra_share, endoorig, exoorig, endobase, exobase, initperiods, simperiods, lastperiods, recompute_final_steady_state, oo_.steady_state, oo_.exo_steady_state); + [extra_completed_share, extra_endo_simul, extra_exo_simul, extra_steady_state, extra_exo_steady_state] = homotopy_loop(extra_share, shareorig, endoorig, exoorig, endobase, exobase, initperiods, simperiods, lastperiods, recompute_final_steady_state, oo_.steady_state, oo_.exo_steady_state); extra_success = (extra_completed_share == extra_share); end extra_simul_time_elapsed = toc(extra_simul_time_counter); @@ -272,7 +294,7 @@ assignin('base', 'Simulated_time_series', ts); oo_.gui.ran_perfect_foresight = oo_.deterministic_simulation.status; -function [completed_share, endo_simul, exo_simul, steady_state, exo_steady_state, iteration, maxerror, solver_iter, per_block_status] = homotopy_loop(max_share, endoorig, exoorig, endobase, exobase, initperiods, simperiods, lastperiods, recompute_final_steady_state, steady_state, exo_steady_state) +function [completed_share, endo_simul, exo_simul, steady_state, exo_steady_state, iteration, maxerror, solver_iter, per_block_status] = homotopy_loop(max_share, shareorig, endoorig, exoorig, endobase, exobase, initperiods, simperiods, lastperiods, recompute_final_steady_state, steady_state, exo_steady_state) % INPUTS % share [double] the share of the shock that we want to simulate % simperiods [vector] period indices of simulation periods (between initial and terminal conditions) @@ -314,16 +336,16 @@ while step > options_.simul.homotopy_min_step_size iter_time_counter = tic; - [steady_success, endo_simul, exo_simul, steady_state, exo_steady_state] = create_scenario(new_share, endoorig, exoorig, endobase, exobase, initperiods, lastperiods, recompute_final_steady_state, endo_simul, steady_state, exo_steady_state); + [steady_success, endo_simul, exo_simul, steady_state, exo_steady_state] = create_scenario(new_share, shareorig, endoorig, exoorig, endobase, exobase, initperiods, lastperiods, recompute_final_steady_state, endo_simul, steady_state, exo_steady_state); if steady_success % At the first iteration, use the initial guess given by - % perfect_foresight_setup or the user (but only if new_share=1, otherwise it + % perfect_foresight_setup or the user (but only if new_share==shareorig, otherwise it % does not make much sense). Afterwards, until a converging iteration has been obtained, % use the rescaled terminal condition (or, if there is no lead, the base % scenario / initial steady state). if completed_share == 0 - if iteration == 1 && new_share == 1 + if iteration == 1 && new_share == shareorig % Nothing to do, at this point endo_simul(:, simperiods) == endoorig(:, simperiods) elseif M_.maximum_lead > 0 endo_simul(:, simperiods) = repmat(endo_simul(:, lastperiods(1)), 1, options_.periods); @@ -409,15 +431,16 @@ end fprintf('\n') -function [steady_success, endo_simul, exo_simul, steady_state, exo_steady_state] = create_scenario(share, endoorig, exoorig, endobase, exobase, initperiods, lastperiods, recompute_final_steady_state, endo_simul, steady_state, exo_steady_state) +function [steady_success, endo_simul, exo_simul, steady_state, exo_steady_state] = create_scenario(share, shareorig, endoorig, exoorig, endobase, exobase, initperiods, lastperiods, recompute_final_steady_state, endo_simul, steady_state, exo_steady_state) % For a given share, comutes the exogenous path and also the initial and % terminal conditions for the endogenous path (but do not modify the initial % guess for endogenous) % % INPUTS % share [double] the share of the shock that we want to simulate -% endoorig [matrix] path of endogenous corresponding to 100% of the shock (only initial and terminal conditions are used) -% exoorig [matrix] path of exogenous corresponding to 100% of the shock +% shareorig [double] the share to which endoorig and exoorig correspond (typically 100%, except for perfect_foresight_with_expectation_errors_solver with homotopy and marginal linearization) +% endoorig [matrix] path of endogenous corresponding to shareorig of the shock (only initial and terminal conditions are used) +% exoorig [matrix] path of exogenous corresponding to shareorig of the shock % endobase [matrix] path of endogenous corresponding to 0% of the shock (only initial and terminal conditions are used) % exobase [matrix] path of exogenous corresponding to 0% of the shock % initperiods [vector] period indices of initial conditions @@ -437,12 +460,12 @@ function [steady_success, endo_simul, exo_simul, steady_state, exo_steady_state] global M_ options_ oo_ % Compute convex combination for the path of exogenous -exo_simul = exoorig*share + exobase*(1-share); +exo_simul = exoorig*share/shareorig + exobase*(1-share); % Compute convex combination for the initial condition % In most cases, the initial condition is a steady state and this does nothing % This is for cases when the initial condition is out of equilibrium -endo_simul(:, initperiods) = share*endoorig(:, initperiods)+(1-share)*endobase(:, initperiods); +endo_simul(:, initperiods) = share/shareorig*endoorig(:, initperiods)+(1-share)*endobase(:, initperiods); % If there is a permanent shock, ensure that the rescaled terminal condition is % a steady state (if the user asked for this recomputation, or if the original @@ -495,7 +518,7 @@ if recompute_final_steady_state options_.markowitz = saved_steady_markowitz; else % The terminal condition is not a steady state, compute a convex combination - endo_simul(:, lastperiods) = share*endoorig(:, lastperiods)+(1-share)*endobase(:, lastperiods); + endo_simul(:, lastperiods) = share/shareorig*endoorig(:, lastperiods)+(1-share)*endobase(:, lastperiods); end diff --git a/matlab/perfect-foresight-models/perfect_foresight_with_expectation_errors_solver.m b/matlab/perfect-foresight-models/perfect_foresight_with_expectation_errors_solver.m index 14ea5a48c..19b3a7211 100644 --- a/matlab/perfect-foresight-models/perfect_foresight_with_expectation_errors_solver.m +++ b/matlab/perfect-foresight-models/perfect_foresight_with_expectation_errors_solver.m @@ -67,7 +67,18 @@ while info_period <= periods options_.periods = sim_length; - perfect_foresight_solver(true); + if info_period > 1 && homotopy_completion_share < 1 && options_.simul.homotopy_marginal_linearization_fallback > 0 + marginal_linearization_previous_raw_sims.sim1.endo_simul = oo_.deterministic_simulation.sim1.endo_simul(:, info_period:end); + marginal_linearization_previous_raw_sims.sim1.exo_simul = oo_.deterministic_simulation.sim1.exo_simul(info_period:end, :); + marginal_linearization_previous_raw_sims.sim1.homotopy_completion_share = oo_.deterministic_simulation.sim1.homotopy_completion_share; + marginal_linearization_previous_raw_sims.sim2.endo_simul = oo_.deterministic_simulation.sim2.endo_simul(:, info_period:end); + marginal_linearization_previous_raw_sims.sim2.exo_simul = oo_.deterministic_simulation.sim2.exo_simul(info_period:end, :); + marginal_linearization_previous_raw_sims.sim2.homotopy_completion_share = oo_.deterministic_simulation.sim2.homotopy_completion_share; + else + marginal_linearization_previous_raw_sims = []; + end + + perfect_foresight_solver(true, marginal_linearization_previous_raw_sims); if ~oo_.deterministic_simulation.status error('perfect_foresight_with_expectation_errors_solver: failed to compute solution for information available at period %d\n', info_period) diff --git a/meson.build b/meson.build index 9ad2a38d4..f8d230fdf 100644 --- a/meson.build +++ b/meson.build @@ -1333,7 +1333,8 @@ mod_and_m_tests = [ 'extra' : [ 'deterministic_simulations/pfwee.csv' ] }, { 'test' : [ 'deterministic_simulations/pfwee_learnt_in.mod' ] }, { 'test' : [ 'deterministic_simulations/pfwee_multiple_shocks.mod' ] }, - { 'test' : [ 'deterministic_simulations/pfwee_homotopy.mod' ] }, + { 'test' : [ 'deterministic_simulations/pfwee_homotopy_linearization.mod' ] }, + { 'test' : [ 'deterministic_simulations/pfwee_homotopy_marginal_linearization.mod' ] }, { 'test' : [ 'lmmcp/rbc.mod' ] }, { 'test' : [ 'lmmcp/sw_lmmcp.mod', 'lmmcp/sw_newton.mod' ], diff --git a/tests/deterministic_simulations/pfwee_homotopy.mod b/tests/deterministic_simulations/pfwee_homotopy_linearization.mod similarity index 95% rename from tests/deterministic_simulations/pfwee_homotopy.mod rename to tests/deterministic_simulations/pfwee_homotopy_linearization.mod index c4175c090..72471830a 100644 --- a/tests/deterministic_simulations/pfwee_homotopy.mod +++ b/tests/deterministic_simulations/pfwee_homotopy_linearization.mod @@ -1,5 +1,5 @@ /* Example that triggers homotopy in perfect foresight simulation with - expectation errors. */ + expectation errors, and tests linearization. */ var Consumption, Capital, LoggedProductivity; diff --git a/tests/deterministic_simulations/pfwee_homotopy_marginal_linearization.mod b/tests/deterministic_simulations/pfwee_homotopy_marginal_linearization.mod new file mode 100644 index 000000000..d884a1029 --- /dev/null +++ b/tests/deterministic_simulations/pfwee_homotopy_marginal_linearization.mod @@ -0,0 +1,40 @@ +/* Example that triggers homotopy in perfect foresight simulation with + expectation errors, and tests marginal linearization. */ + +var Consumption, Capital, LoggedProductivity; + +varexo LoggedProductivityInnovation; + +parameters beta, alpha, delta, rho; + +beta = .985; +alpha = 1/3; +delta = alpha/10; +rho = .9; + +model; + 1/Consumption = beta/Consumption(1)*(alpha*exp(LoggedProductivity(1))*Capital^(alpha-1)+1-delta); + Capital = exp(LoggedProductivity)*Capital(-1)^alpha+(1-delta)*Capital(-1)-Consumption; + LoggedProductivity = rho*LoggedProductivity(-1)+LoggedProductivityInnovation; +end; + +initval; + LoggedProductivityInnovation = 0; +end; + +steady; + +endval; + LoggedProductivityInnovation = 0.6; +end; + +endval(learnt_in = 5); + LoggedProductivityInnovation = 1; +end; + +perfect_foresight_with_expectation_errors_setup(periods=200); +perfect_foresight_with_expectation_errors_solver(homotopy_max_completion_share = 0.8, homotopy_marginal_linearization_fallback, steady_solve_algo = 13); + +if ~oo_.deterministic_simulation.status + error('Perfect foresight simulation failed') +end