From 24cc67e585051f8ec6f3351b930397f9ea333ce6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Fri, 26 Apr 2019 15:36:33 +0200 Subject: [PATCH] Ensure that all perfect foresight solvers work with periods=1. See #1205 and #1176. --- .../linear_perfect_foresight_problem.m | 29 ++++---- .../perfect_foresight_mcp_problem.m | 13 ++-- .../perfect_foresight_problem.m | 72 +++++++++---------- .../perfect_foresight_solver.m | 27 ++++--- .../perfect_foresight_solver_core.m | 41 ++++++----- .../private/initialize_stacked_problem.m | 20 ++++-- .../private/simulation_core.m | 31 ++++---- matlab/perfect-foresight-models/sim1.m | 2 +- matlab/perfect-foresight-models/sim1_linear.m | 2 +- .../solve_stacked_linear_problem.m | 4 +- .../solve_stacked_problem.m | 11 ++- 11 files changed, 142 insertions(+), 110 deletions(-) diff --git a/matlab/perfect-foresight-models/linear_perfect_foresight_problem.m b/matlab/perfect-foresight-models/linear_perfect_foresight_problem.m index 3a04ea8e1..fd76cd89e 100644 --- a/matlab/perfect-foresight-models/linear_perfect_foresight_problem.m +++ b/matlab/perfect-foresight-models/linear_perfect_foresight_problem.m @@ -1,25 +1,22 @@ 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,nnzJ,jendo,jexog) -% function [residuals,JJacobian] = perfect_foresight_problem(x, model_dynamic, Y0, YT,exo_simul, -% params, steady_state, maximum_lag, periods, ny, i_cols,i_cols_J1, i_cols_1, -% i_cols_T, i_cols_j, nnzA) -% computes the residuals and th Jacobian matrix -% for a perfect foresight problem over T periods. + 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, nnzJ, jendo, jexog) + +% Computes the residuals and the Jacobian matrix for a linear perfect foresight problem over T periods. % % INPUTS -% ... +% ... +% % OUTPUTS -% ... +% ... +% % ALGORITHM -% ... +% ... % % SPECIAL REQUIREMENTS % None. -% Copyright (C) 2015-2017 Dynare Team +% Copyright (C) 2015-2019 Dynare Team % % This file is part of Dynare. % @@ -44,7 +41,7 @@ residuals = zeros(T*ny,1); z = zeros(columns(dynamicjacobian), 1); if nargout == 2 - JJacobian = sparse([],[],[],T*ny,T*ny,T*nnzJ); + JJacobian = spalloc(T*ny, T*ny, T*nnzJ); end i_rows = 1:ny; @@ -55,7 +52,9 @@ for it = maximum_lag+(1:T) z(jexog) = transpose(exo_simul(it,:)); residuals(i_rows) = dynamicjacobian*z; if nargout == 2 - if it == maximum_lag+1 + if T==1 && it==maximum_lag+1 + JJacobian(i_rows, i_cols_J0) = dynamicjacobian(:,i_cols_0); + elseif it == maximum_lag+1 JJacobian(i_rows,i_cols_J1) = dynamicjacobian(:,i_cols_1); elseif it == maximum_lag+T JJacobian(i_rows,i_cols_J(i_cols_T)) = dynamicjacobian(:,i_cols_T); diff --git a/matlab/perfect-foresight-models/perfect_foresight_mcp_problem.m b/matlab/perfect-foresight-models/perfect_foresight_mcp_problem.m index 000fff50f..89d9a7b82 100644 --- a/matlab/perfect-foresight-models/perfect_foresight_mcp_problem.m +++ b/matlab/perfect-foresight-models/perfect_foresight_mcp_problem.m @@ -2,7 +2,7 @@ function [residuals,JJacobian] = perfect_foresight_mcp_problem(y, dynamic_functi exo_simul, params, steady_state, ... maximum_lag, T, ny, i_cols, ... i_cols_J1, i_cols_1, i_cols_T, ... - i_cols_j,nnzJ,eq_index) + i_cols_j, i_cols_0,i_cols_J0, nnzJ,eq_index) % function [residuals,JJacobian] = perfect_foresight_mcp_problem(y, dynamic_function, Y0, YT, ... % exo_simul, params, steady_state, ... % maximum_lag, T, ny, i_cols, ... @@ -80,10 +80,12 @@ for it = maximum_lag+(1:T) steady_state,it); residuals(i_rows) = res(eq_index); elseif nargout == 2 - [res,jacobian] = dynamic_function(YY(i_cols),exo_simul, params, ... - steady_state,it); + [res,jacobian] = dynamic_function(YY(i_cols),exo_simul, params, steady_state,it); residuals(i_rows) = res(eq_index); - if it == maximum_lag+1 + 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(eq_index,i_cols_1)); iJacobian{1} = [offset+rows, i_cols_J1(cols), vals]; elseif it == maximum_lag+T @@ -103,6 +105,5 @@ end if nargout == 2 iJacobian = cat(1,iJacobian{:}); - JJacobian = sparse(iJacobian(:,1),iJacobian(:,2),iJacobian(:,3),T* ... - ny,T*ny); + JJacobian = sparse(iJacobian(:,1),iJacobian(:,2),iJacobian(:,3),T*ny,T*ny); end \ No newline at end of file diff --git a/matlab/perfect-foresight-models/perfect_foresight_problem.m b/matlab/perfect-foresight-models/perfect_foresight_problem.m index c6c6f77a7..e71c25248 100644 --- a/matlab/perfect-foresight-models/perfect_foresight_problem.m +++ b/matlab/perfect-foresight-models/perfect_foresight_problem.m @@ -2,47 +2,48 @@ function [residuals,JJacobian] = perfect_foresight_problem(y, dynamic_function, exo_simul, params, steady_state, ... maximum_lag, T, ny, i_cols, ... i_cols_J1, i_cols_1, i_cols_T, ... - i_cols_j,nnzJ) -% 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,nnzJ) -% computes the residuals and the Jacobian matrix for a perfect foresight problem over T periods. + i_cols_j, i_cols_0, i_cols_J0, nnzJ) + +% 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) +% - 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 +% - 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 +% - 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 +% - 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 +% - 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 +% - i_cols_j [double] indices of contemporaneous variables in M.lead_lag_incidence % in dynamic Jacobian (relevant in intermediate periods) -% nnzJ [scalar] number of non-zero elements in Jacobian +% - 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) +% - nnzJ [scalar] number of non-zero elements in Jacobian +% % OUTPUTS -% residuals [double] (N*T)*1 array, residuals of the stacked problem -% JJacobian [double] (N*T)*(N*T) array, Jacobian of the stacked problem +% - 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 +% None % % SPECIAL REQUIREMENTS -% None. +% None. -% Copyright (C) 1996-2017 Dynare Team +% Copyright (C) 1996-2019 Dynare Team % % This file is part of Dynare. % @@ -73,12 +74,13 @@ 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); + 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 it == maximum_lag+1 + [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 @@ -91,13 +93,11 @@ for it = maximum_lag+(1:T) 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); + JJacobian = sparse(iJacobian(:,1), iJacobian(:,2), iJacobian(:,3), T*ny, T*ny); end \ No newline at end of file diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver.m b/matlab/perfect-foresight-models/perfect_foresight_solver.m index 17f607ca0..c7fab2a6f 100644 --- a/matlab/perfect-foresight-models/perfect_foresight_solver.m +++ b/matlab/perfect-foresight-models/perfect_foresight_solver.m @@ -37,6 +37,8 @@ if isempty(options_.scalv) || options_.scalv == 0 options_.scalv = oo_.steady_state; end +periods = options_.periods; + options_.scalv= 1; if options_.debug @@ -56,7 +58,7 @@ if options_.debug end initperiods = 1:M_.maximum_lag; -lastperiods = (M_.maximum_lag+options_.periods+1):(M_.maximum_lag+options_.periods+M_.maximum_lead); +lastperiods = (M_.maximum_lag+periods+1):(M_.maximum_lag+periods+M_.maximum_lead); oo_ = perfect_foresight_solver_core(M_,options_,oo_); @@ -91,8 +93,8 @@ if ~oo_.deterministic_simulation.status && ~options_.no_homotopy options_.verbosity = 0; % Set initial paths for the endogenous and exogenous variables. - endoinit = repmat(oo_.steady_state, 1,M_.maximum_lag+options_.periods+M_.maximum_lead); - exoinit = repmat(oo_.exo_steady_state',M_.maximum_lag+options_.periods+M_.maximum_lead,1); + endoinit = repmat(oo_.steady_state, 1,M_.maximum_lag+periods+M_.maximum_lead); + exoinit = repmat(oo_.exo_steady_state',M_.maximum_lag+periods+M_.maximum_lead,1); % Copy the current paths for the exogenous and endogenous variables. exosim = oo_.exo_simul; @@ -131,7 +133,7 @@ if ~oo_.deterministic_simulation.status && ~options_.no_homotopy if isequal(iteration, 1) % First iteration, same initial guess as in the first call to perfect_foresight_solver_core routine. - oo_.endo_simul(:,M_.maximum_lag+1:end-M_.maximum_lead) = endoinit(:,1:options_.periods); + oo_.endo_simul(:,M_.maximum_lag+1:end-M_.maximum_lead) = endoinit(:,1:periods); elseif path_with_nans || path_with_cplx % If solver failed with NaNs or complex number, use previous solution as an initial guess. oo_.endo_simul(:,M_.maximum_lag+1:end-M_.maximum_lead) = saved_endo_simul(:,1+M_.maximum_lag:end-M_.maximum_lead); @@ -174,19 +176,26 @@ if ~oo_.deterministic_simulation.status && ~options_.no_homotopy end -if ~isreal(oo_.endo_simul(:)) %can only happen without bytecode +if ~isreal(oo_.endo_simul(:)) % can only happen without bytecode y0 = real(oo_.endo_simul(:,1)); - yT = real(oo_.endo_simul(:,options_.periods+2)); - yy = real(oo_.endo_simul(:,2:options_.periods+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,:)'); + else + i_cols_0 = []; + i_cols_J0 = []; + end residuals = perfect_foresight_problem(yy(:),str2func([M_.fname '.dynamic']), y0, yT, ... oo_.exo_simul,M_.params,oo_.steady_state, ... - M_.maximum_lag,options_.periods,M_.endo_nbr,i_cols, ... - i_cols_J1, i_cols_1, i_cols_T, i_cols_j, ... + 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, ... M_.NNZDerivatives(1)); if max(abs(residuals))< options_.dynatol.f oo_.deterministic_simulation.status = 1; diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m index 1f55a2de4..f0a73f580 100644 --- a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m +++ b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m @@ -1,5 +1,5 @@ function [oo_, maxerror] = perfect_foresight_solver_core(M_, options_, oo_) -%function [oo_, maxerror] = perfect_foresight_solver_core(M_, options_, oo_) + % Core function calling solvers for perfect foresight model % % INPUTS @@ -11,7 +11,7 @@ function [oo_, maxerror] = perfect_foresight_solver_core(M_, options_, oo_) % - oo_ [struct] contains results % - maxerror [double] contains the maximum absolute error -% Copyright (C) 2015-2017 Dynare Team +% Copyright (C) 2015-2019 Dynare Team % % This file is part of Dynare. % @@ -33,18 +33,20 @@ if options_.lmmcp.status options_.solve_algo = 10; end +periods = options_.periods; + if options_.linear_approximation && ~(isequal(options_.stack_solve_algo,0) || isequal(options_.stack_solve_algo,7)) - error('perfect_foresight_solver: Option linear_approximation is only available with option stack_solve_algo equal to 0.') + error('perfect_foresight_solver: Option linear_approximation is only available with option stack_solve_algo equal to 0 or 7.') end -if options_.linear && isequal(options_.stack_solve_algo,0) - options_.linear_approximation = 1; +if options_.linear && (isequal(options_.stack_solve_algo, 0) || isequal(options_.stack_solve_algo, 7)) + options_.linear_approximation = true; end if options_.block if options_.bytecode try - [info, tmp] = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state,1,options_.periods+2), options_.periods); + [info, tmp] = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state,1, periods+2), periods); catch info = 1; end @@ -63,7 +65,7 @@ if options_.block else if options_.bytecode try - [info, tmp] = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state,1,options_.periods+2), options_.periods); + [info, tmp] = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state, 1, periods+2), periods); catch info = 1; end @@ -119,14 +121,19 @@ end if nargout>1 y0 = oo_.endo_simul(:,1); - yT = oo_.endo_simul(:,options_.periods+2); - yy = oo_.endo_simul(:,2:options_.periods+1); - if ~exist('illi') - 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,:)'); + 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; @@ -136,8 +143,8 @@ if nargout>1 else residuals = perfect_foresight_problem(yy(:),str2func([M_.fname '.dynamic']), y0, yT, ... oo_.exo_simul,M_.params,oo_.steady_state, ... - M_.maximum_lag,options_.periods,M_.endo_nbr,i_cols, ... - i_cols_J1, i_cols_1, i_cols_T, i_cols_j, ... + 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, ... M_.NNZDerivatives(1)); end maxerror = max(max(abs(residuals))); diff --git a/matlab/perfect-foresight-models/private/initialize_stacked_problem.m b/matlab/perfect-foresight-models/private/initialize_stacked_problem.m index f019805dc..55290fc9f 100644 --- a/matlab/perfect-foresight-models/private/initialize_stacked_problem.m +++ b/matlab/perfect-foresight-models/private/initialize_stacked_problem.m @@ -1,7 +1,6 @@ -function [options, y0, yT, z, i_cols, i_cols_J1, i_cols_T, i_cols_j, i_cols_1, ... - dynamicmodel] = initialize_stacked_problem(endogenousvariables, options, M, steadystate_y) -% function [options, y0, yT, z, i_cols, i_cols_J1, i_cols_T, i_cols_j, i_cols_1, ... -% dynamicmodel] = initialize_stacked_problem(endogenousvariables, options, M, steadystate_y) +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 @@ -9,6 +8,7 @@ function [options, y0, yT, z, i_cols, i_cols_J1, i_cols_T, i_cols_j, i_cols_1, . % - 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 @@ -25,9 +25,12 @@ function [options, y0, yT, z, i_cols, i_cols_J1, i_cols_T, i_cols_j, i_cols_1, . % 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 (C) 2015-2017 Dynare Team +% Copyright (C) 2015-2019 Dynare Team % % This file is part of Dynare. % @@ -75,4 +78,11 @@ illi = M.lead_lag_incidence'; 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 dynamicmodel = str2func([M.fname,'.dynamic']); \ No newline at end of file diff --git a/matlab/perfect-foresight-models/private/simulation_core.m b/matlab/perfect-foresight-models/private/simulation_core.m index aa66ef2de..6e7bae6a0 100644 --- a/matlab/perfect-foresight-models/private/simulation_core.m +++ b/matlab/perfect-foresight-models/private/simulation_core.m @@ -1,7 +1,7 @@ function [oo_, maxerror] = simulation_core(options_, M_, oo_) %function [oo_, maxerror] = simulation_core(options_, M_, oo_) -% Copyright (C) 2015-2017 Dynare Team +% Copyright (C) 2015-2019 Dynare Team % % This file is part of Dynare. % @@ -26,10 +26,12 @@ if options_.linear && isequal(options_.stack_solve_algo,0) options_.linear_approximation = 1; end +periods = options_.periods; + if options_.block if options_.bytecode try - [info, tmp] = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state,1,options_.periods+2), options_.periods); + [info, tmp] = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state, 1, periods+2), periods); catch info = 0; end @@ -48,7 +50,7 @@ if options_.block else if options_.bytecode try - [info, tmp] = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state,1,options_.periods+2), options_.periods); + [info, tmp] = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state, 1, periods+2), periods); catch info = 0; end @@ -95,14 +97,19 @@ end if nargout>1 y0 = oo_.endo_simul(:,1); - yT = oo_.endo_simul(:,options_.periods+2); - yy = oo_.endo_simul(:,2:options_.periods+1); - if ~exist('illi') - 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,:)'); + 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; @@ -113,7 +120,7 @@ if nargout>1 residuals = perfect_foresight_problem(yy(:),str2func([M_.fname '.dynamic']), y0, yT, ... oo_.exo_simul,M_.params,oo_.steady_state, ... M_.maximum_lag,options_.periods,M_.endo_nbr,i_cols, ... - i_cols_J1, i_cols_1, i_cols_T, i_cols_j, ... + i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0, ... M_.NNZDerivatives(1)); end maxerror = max(max(abs(residuals))); diff --git a/matlab/perfect-foresight-models/sim1.m b/matlab/perfect-foresight-models/sim1.m index f43d6df0d..66e698c7c 100644 --- a/matlab/perfect-foresight-models/sim1.m +++ b/matlab/perfect-foresight-models/sim1.m @@ -101,7 +101,7 @@ for iter = 1:options.simul.maxit m = 0; for it = (maximum_lag+1):(maximum_lag+periods) [d1,jacobian] = model_dynamic(Y(i_cols), exogenousvariables, params, steadystate,it); - if it == maximum_lag+periods && it == maximum_lag+1 + 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 diff --git a/matlab/perfect-foresight-models/sim1_linear.m b/matlab/perfect-foresight-models/sim1_linear.m index 3fa40bf41..2a68ac0b5 100644 --- a/matlab/perfect-foresight-models/sim1_linear.m +++ b/matlab/perfect-foresight-models/sim1_linear.m @@ -142,7 +142,7 @@ i_cols_A = ipcn; i_cols = ipcn+(maximum_lag-1)*ny; m = 0; for it = (maximum_lag+1):(maximum_lag+periods) - if isequal(it, maximum_lag+periods) && isequal(it, maximum_lag+1) + if periods==1 && isequal(it, maximum_lag+1) nv = length(v0); iA(iv0+m,:) = [i_rows(r0),ic(c0),v0]; elseif isequal(it, maximum_lag+periods) diff --git a/matlab/perfect-foresight-models/solve_stacked_linear_problem.m b/matlab/perfect-foresight-models/solve_stacked_linear_problem.m index 5464c9e06..cd52ac444 100644 --- a/matlab/perfect-foresight-models/solve_stacked_linear_problem.m +++ b/matlab/perfect-foresight-models/solve_stacked_linear_problem.m @@ -17,7 +17,7 @@ function [endogenousvariables, info] = solve_stacked_linear_problem(endogenousva % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -[options, y0, yT, z, i_cols, i_cols_J1, i_cols_T, i_cols_j, i_cols_1, dynamicmodel] = ... +[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,:)'); @@ -45,7 +45,7 @@ x = bsxfun(@minus, exogenousvariables, steadystate_x'); 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_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0, ... M.NNZDerivatives(1), jendo, jexog); if all(imag(y)<.1*options.dynatol.x) diff --git a/matlab/perfect-foresight-models/solve_stacked_problem.m b/matlab/perfect-foresight-models/solve_stacked_problem.m index 3a1e2e715..d673a1973 100644 --- a/matlab/perfect-foresight-models/solve_stacked_problem.m +++ b/matlab/perfect-foresight-models/solve_stacked_problem.m @@ -1,5 +1,5 @@ function [endogenousvariables, info] = solve_stacked_problem(endogenousvariables, exogenousvariables, steadystate, M, options) -% [endogenousvariables, info] = solve_stacked_problem(endogenousvariables, exogenousvariables, steadystate, M, options) + % Solves the perfect foresight model using dynare_solve % % INPUTS @@ -13,7 +13,7 @@ function [endogenousvariables, info] = solve_stacked_problem(endogenousvariables % - endogenousvariables [double] N*T array, paths for the endogenous variables (solution of the perfect foresight model). % - info [struct] contains informations about the results. -% Copyright (C) 2015-2017 Dynare Team +% Copyright (C) 2015-2019 Dynare Team % % This file is part of Dynare. % @@ -30,7 +30,7 @@ function [endogenousvariables, info] = solve_stacked_problem(endogenousvariables % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -[options, y0, yT, z, i_cols, i_cols_J1, i_cols_T, i_cols_j, i_cols_1, dynamicmodel] = ... +[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); if (options.solve_algo == 10 || options.solve_algo == 11)% mixed complementarity problem @@ -50,14 +50,14 @@ if (options.solve_algo == 10 || options.solve_algo == 11)% mixed complementarity dynamicmodel, y0, yT, ... exogenousvariables, M.params, steadystate, ... M.maximum_lag, options.periods, M.endo_nbr, i_cols, ... - i_cols_J1, i_cols_1, i_cols_T, i_cols_j, ... + i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0, ... M.NNZDerivatives(1),eq_index); else [y, check] = dynare_solve(@perfect_foresight_problem,z(:),options, ... dynamicmodel, y0, yT, ... exogenousvariables, M.params, steadystate, ... M.maximum_lag, options.periods, M.endo_nbr, i_cols, ... - i_cols_J1, i_cols_1, i_cols_T, i_cols_j, ... + i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0, ... M.NNZDerivatives(1)); end @@ -74,6 +74,5 @@ endogenousvariables(:, M.maximum_lag+(1:options.periods)) = reshape(y, M.endo_nb if check info.status = false; else - info.status = true; end \ No newline at end of file