Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • giovanma/dynare
  • giorgiomas/dynare
  • Vermandel/dynare
  • Dynare/dynare
  • normann/dynare
  • MichelJuillard/dynare
  • wmutschl/dynare
  • FerhatMihoubi/dynare
  • sebastien/dynare
  • lnsongxf/dynare
  • rattoma/dynare
  • CIMERS/dynare
  • FredericKarame/dynare
  • SumuduK/dynare
  • MinjeJeon/dynare
  • camilomrch/dynare
  • DoraK/dynare
  • avtishin/dynare
  • selma/dynare
  • claudio_olguin/dynare
  • jeffjiang07/dynare
  • EthanSystem/dynare
  • stepan-a/dynare
  • wjgatt/dynare
  • JohannesPfeifer/dynare
  • gboehl/dynare
26 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 698 additions and 456 deletions
matlab/discretionary_policy/discretionary_policy.m
View file @ e9d4a40a
......@@ -35,6 +35,9 @@ origorder = options_.order;
options_.order = 1;
[info, oo_, options_, M_] = stoch_simul(M_, options_, oo_, var_list);
if info(1)
return;
end
oo_.steady_state = oo_.dr.ys;
if ~options_.noprint
......
matlab/discretionary_policy/discretionary_policy_1.m
View file @ e9d4a40a
......@@ -42,6 +42,7 @@ if options_.steadystate_flag
[ys,M_.params,info] = evaluate_steady_state_file(endo_steady_state,[exo_steady_state; exo_det_steady_state],M_, ...
options_,~options_.steadystate.nocheck);
if info(1)
params=M_.params;
return;
end
else
......@@ -113,5 +114,7 @@ T=H(dr.order_var,dr.order_var);
dr.ghu=G(dr.order_var,:);
if M_.maximum_endo_lag
Selection=M_.lead_lag_incidence(1,dr.order_var)>0;%select state variables
else
Selection=[];
end
dr.ghx=T(:,Selection);
matlab/dsample.m 0 → 100644
View file @ e9d4a40a
function dsample(s1,s2)
% function dsample(s1,s2)
% This optional command permits to reduce the number of periods considered in following output commands.
% If only one argument is provided, output is from period 1 to the period specified in the DSAMPLE command.
% If two arguments are present output is done for the interval between the two periods.
% DSAMPLE without arguments reset the sample to the one specified by PERIODS
%
% INPUTS
% s1: first period
% s2: last period
%
% OUTPUTS
% none
%
% SPECIAL REQUIREMENTS
% none
% Note: this is a standalone command documented in the manual
% Copyright © 2001-2025 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/>.
global options_
options_.smpl = zeros(2,1) ;
if nargin == 0
options_.smpl(1) = 1 ;
options_.smpl(2) = options_.periods ;
elseif nargin == 1
if s1 > options_.periods
error('DSAMPLE: argument greater than number of periods');
end
options_.smpl(1) = 1 ;
options_.smpl(2) = s1 ;
else
if s1 > options_.periods || s2 > options_.periods
error('DSAMPLE: one of the arguments is greater than number of periods');
end
options_.smpl(1) = s1 ;
options_.smpl(2) = s2 ;
end
% 02/23/01 MJ added error checking
\ No newline at end of file
dseries @ 657f8b8e
Compare a456425a...657f8b8e
Subproject commit a456425a6ff65c699f5937c0ea775de1a3348405
Subproject commit 657f8b8eaa4968c6f28434822d0f063f9c955935
matlab/dynare_config.m
View file @ e9d4a40a
......@@ -53,6 +53,7 @@ p = {'/../contrib/ms-sbvar/TZcode/MatlabFiles/' ; ...
'/ep/' ; ...
'/estimation/'; ...
'/estimation/smc/'; ...
'/estimation/online/'; ...
'/estimation/resampler/'; ...
'/kalman/' ; ...
'/kalman/likelihood' ; ...
......
matlab/ep/ep_problem_0.m 0 → 100644
View file @ e9d4a40a
function [res, A, info] = ep_problem_0(y, x, pfm)
% Evaluate the residuals and stacked jacobian of a stochastic perfect
% foresight, considering sequences of future innovations in a perfect n-ary tree.
%
% INPUTS:
% - y [double] m×1 vector (endogenous variables in all periods and future worlds).
% - x [double] q×1 vector of exogenous variables.
% - pfm [struct] Definition of the perfect foresight model to be solved.
%
% OUTPUTS:
% - res [double] m×1 vector, residuals of the stacked equations.
% - A [double] m×m sparse matrix, jacobian of the stacked equations.
% - info [logical] scalar
%
% REMARKS:
% [1] The structure pfm holds the given initial condition for the states (pfm.y0) and the terminal condition
info = false;
params = pfm.params;
steady_state = pfm.steady_state;
ny = pfm.ny;
periods = pfm.periods;
dynamic_model = pfm.dynamic_model;
lead_lag_incidence = pfm.lead_lag_incidence;
i_cols_1 = pfm.i_cols_1;
i_cols_j = pfm.i_cols_j;
i_cols_T = pfm.i_cols_T;
order = pfm.order;
hybrid_order = pfm.hybrid_order;
h_correction = pfm.h_correction;
nodes = pfm.nodes;
weights = pfm.weights;
nnodes = pfm.nnodes;
i_cols_p = pfm.i_cols_p;
i_cols_s = pfm.i_cols_s;
i_cols_f = pfm.i_cols_f;
i_rows = pfm.i_rows;
i_cols_Ap = pfm.i_cols_Ap;
i_cols_As = pfm.i_cols_As;
i_cols_Af = pfm.i_cols_Af;
i_hc = pfm.i_hc;
Y = pfm.Y;
Y(pfm.i_upd_y) = y;
A1 = pfm.A1;
res = pfm.res;
for i = 1:order+1
i_w_p = 1;
for j = 1:nnodes^(i-1)
innovation = x;
if i > 1
innovation(i+1,:) = nodes(mod(j-1,nnodes)+1,:);
end
if i <= order
for k=1:nnodes
if hybrid_order && i==order
z = [Y(i_cols_p,i_w_p);
Y(i_cols_s,j);
Y(i_cols_f,(j-1)*nnodes+k)+h_correction(i_hc)];
else
z = [Y(i_cols_p,i_w_p);
Y(i_cols_s,j);
Y(i_cols_f,(j-1)*nnodes+k)];
end
[d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1);
if i == 1
% in first period we don't keep track of
% predetermined variables
i_cols_A = [i_cols_As - ny; i_cols_Af];
A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_1);
else
i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af];
A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_j);
end
res(:,i,j) = res(:,i,j)+weights(k)*d1;
i_cols_Af = i_cols_Af + ny;
end
else
z = [Y(i_cols_p,i_w_p);
Y(i_cols_s,j);
Y(i_cols_f,j)];
[d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1);
if i == 1
% in first period we don't keep track of
% predetermined variables
i_cols_A = [i_cols_As - ny; i_cols_Af];
A1(i_rows,i_cols_A) = jacobian(:,i_cols_1);
else
i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af];
A1(i_rows,i_cols_A) = jacobian(:,i_cols_j);
end
res(:,i,j) = d1;
i_cols_Af = i_cols_Af + ny;
end
i_rows = i_rows + ny;
if mod(j,nnodes) == 0
i_w_p = i_w_p + 1;
end
if i > 1
if mod(j,nnodes) == 0
i_cols_Ap = i_cols_Ap + ny;
end
i_cols_As = i_cols_As + ny;
end
end
i_cols_p = i_cols_p + ny;
i_cols_s = i_cols_s + ny;
i_cols_f = i_cols_f + ny;
end
nzA = cell(periods,pfm.world_nbr);
for j=1:pfm.world_nbr
i_rows_y = find(lead_lag_incidence')+(order+1)*ny;
offset_c = ny*(sum(nnodes.^(0:order-1),2)+j-1);
offset_r = (j-1)*ny;
for i=order+2:periods
[d1,jacobian] = dynamic_model(Y(i_rows_y,j), x, params, steady_state, i+1);
if i == periods
[ir,ic,v] = find(jacobian(:,i_cols_T));
else
[ir,ic,v] = find(jacobian(:,i_cols_j));
end
nzA{i,j} = [offset_r+ir,offset_c+pfm.icA(ic), v]';
res(:,i,j) = d1;
i_rows_y = i_rows_y + ny;
offset_c = offset_c + pfm.world_nbr*ny;
offset_r = offset_r + pfm.world_nbr*ny;
end
end
A2 = [nzA{:}]';
A = [A1; sparse(A2(:,1),A2(:,2),A2(:,3),ny*(periods-order-1)*pfm.world_nbr,pfm.dimension)];
res = res(pfm.i_upd_r);
matlab/ep/ep_problem_2.m → matlab/ep/ep_problem_1.m
View file @ e9d4a40a
function [res,A,info] = ep_problem_2(y,x,pm)
function [res, A, info] = ep_problem_1(y, x, pfm)
info = 0;
res = [];
% Evaluate the residuals and stacked jacobian of a stochastic perfect
% foresight, considering sequences of future innovations in a sparse tree.
%
% INPUTS:
% - y [double] m×1 vector (endogenous variables in all periods and future worlds).
% - x [double] q×1 vector of exogenous variables.
% - pfm [struct] Definition of the perfect foresight model to be solved.
%
% OUTPUTS:
% - res [double] m×1 vector, residuals of the stacked equations.
% - A [double] m×m sparse matrix, jacobian of the stacked equations.
% - info [logical] scalar
%
% REMARKS:
% [1] The structure pfm holds the given initial condition for the states (pfm.y0) and the terminal condition
info = false;
A = [];
dynamic_model = pm.dynamic_model;
ny = pm.ny;
params = pm.params;
steady_state = pm.steady_state;
order = pm.order;
nodes = pm.nodes;
nnodes = pm.nnodes;
weights = pm.weights;
h_correction = pm.h_correction;
dimension = pm.dimension;
world_nbr = pm.world_nbr;
nnzA = pm.nnzA;
periods = pm.periods;
i_rows = pm.i_rows';
i_cols = pm.i_cols;
nyp = pm.nyp;
nyf = pm.nyf;
hybrid_order = pm.hybrid_order;
i_cols_1 = pm.i_cols_1;
i_cols_j = pm.i_cols_j;
icA = pm.icA;
i_cols_T = pm.i_cols_T;
eq_index = pm.eq_index;
dynamic_model = pfm.dynamic_model;
ny = pfm.ny;
params = pfm.params;
steady_state = pfm.steady_state;
order = pfm.stochastic_order;
nodes = pfm.nodes;
nnodes = pfm.nnodes;
weights = pfm.weights;
h_correction = pfm.h_correction;
dimension = pfm.dimension;
world_nbr = pfm.world_nbr;
nnzA = pfm.nnzA;
periods = pfm.periods;
i_rows = pfm.i_rows';
i_cols = pfm.i_cols;
nyp = pfm.nyp;
nyf = pfm.nyf;
hybrid_order = pfm.hybrid_order;
i_cols_1 = pfm.i_cols_1;
i_cols_j = pfm.i_cols_j;
icA = pfm.icA;
i_cols_T = pfm.i_cols_T;
eq_index = pfm.eq_index;
i_cols_p = i_cols(1:nyp);
i_cols_s = i_cols(nyp+(1:ny));
i_cols_f = i_cols(nyp+ny+(1:nyf));
i_cols_A = i_cols;
i_cols_Ap0 = i_cols_p;
i_cols_As = i_cols_s;
i_cols_Af0 = i_cols_f - ny;
......@@ -40,9 +56,10 @@ i_hc = i_cols_f - 2*ny;
nzA = cell(periods,world_nbr);
res = zeros(ny,periods,world_nbr);
Y = zeros(ny*(periods+2),world_nbr);
Y(1:ny,1) = pm.y0;
Y(1:ny,1) = pfm.y0;
Y(end-ny+1:end,:) = repmat(steady_state,1,world_nbr);
Y(pm.i_upd_y) = y;
Y(pfm.i_upd_y) = y;
offset_r0 = 0;
for i = 1:order+1
i_w_p = 1;
......@@ -78,10 +95,10 @@ for i = 1:order+1
else
k1 = (nnodes-1)*(i-1)+k;
end
if hybrid_order == 2 && (k > 1 || i == order)
if hybrid_order && (k > 1 || i == order)
z = [Y(i_cols_p,1);
Y(i_cols_s,1);
Y(i_cols_f,k1)+h_correction(i_hc)];
Y(i_cols_f,k1)+h_correction(i_hc)];
else
z = [Y(i_cols_p,1);
Y(i_cols_s,1);
......@@ -192,4 +209,4 @@ if nargout > 1
iA = [nzA{:}]';
A = sparse(iA(:,1),iA(:,2),iA(:,3),dimension,dimension);
end
res = res(pm.i_upd_r);
\ No newline at end of file
res = res(pfm.i_upd_r);
matlab/ep/euler_equation_error.m
View file @ e9d4a40a
function e = euler_equation_error(y0,x,innovations,M_,options_,oo_,pfm,nodes,weights)
% e = euler_equation_error(y0,x,innovations,M_,options_,oo_,pfm,nodes,weights)
% Copyright © 2016-2023 Dynare Team
% Copyright © 2016-2025 Dynare Team
%
% This file is part of Dynare.
%
......@@ -23,7 +23,7 @@ ep = options_.ep;
y1 = extended_path_core(ep.periods, ...
M_.endo_nbr, M_.exo_nbr, ...
innovations.positive_var_indx, ...
x, ep.init, y0, oo_.steady_state, ...
x, ep.use_first_order_solution_as_initial_guess, y0, oo_.steady_state, ...
0, ...
ep.stochastic.order, M_, ...
pfm, ep.stochastic.algo, ...
......@@ -38,7 +38,7 @@ for i=1:length(nodes)
x2 = x1;
x2(2,:) = x2(2,:) + nodes(i,:);
y2 = extended_path_core(ep.periods, M_.endo_nbr, M_.exo_nbr, ...
innovations.positive_var_indx, x2, ep.init, ...
innovations.positive_var_indx, x2, ep.use_first_order_solution_as_initial_guess, ...
y1, oo_.steady_state, 0, ...
ep.stochastic.order, M_, pfm, ep.stochastic.algo, ...
ep.solve_algo, ep.stack_solve_algo, options_.lmmcp, ...
......
matlab/ep/extended_path.m
View file @ e9d4a40a
......@@ -19,7 +19,7 @@ function [ts,oo_] = extended_path(initialconditions, samplesize, exogenousvariab
%
% SPECIAL REQUIREMENTS
% Copyright © 2009-2023 Dynare Team
% Copyright © 2009-2025 Dynare Team
%
% This file is part of Dynare.
%
......@@ -36,18 +36,22 @@ function [ts,oo_] = extended_path(initialconditions, samplesize, exogenousvariab
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
[initialconditions, innovations, pfm, ep, verbosity, options_, oo_] = ...
if ~isempty(M_.perfect_foresight_controlled_paths)
error('extended_path command is not compatible with perfect_foresight_controlled_paths block')
end
[initialconditions, innovations, pfm, options_, oo_] = ...
extended_path_initialization(initialconditions, samplesize, exogenousvariables, options_, M_, oo_);
[shocks, spfm_exo_simul, innovations, oo_] = extended_path_shocks(innovations, ep, exogenousvariables, samplesize,M_,options_,oo_);
[shocks, spfm_exo_simul, innovations, oo_] = extended_path_shocks(innovations, exogenousvariables, samplesize, M_, options_, oo_);
% Initialize the matrix for the paths of the endogenous variables.
endogenous_variables_paths = NaN(M_.endo_nbr,samplesize+1);
endogenous_variables_paths = NaN(M_.endo_nbr, samplesize+1);
endogenous_variables_paths(:,1) = initialconditions;
% Set waitbar (graphic or text mode)
hh_fig = dyn_waitbar(0,'Please wait. Extended Path simulations...');
set(hh_fig,'Name','EP simulations.');
set(hh_fig,'Name', 'EP simulations.' );
% Initialize while-loop index.
t = 1;
......@@ -67,14 +71,27 @@ while (t <= samplesize)
else
initialguess = [];
end
[endogenous_variables_paths(:,t), info_convergence, endogenousvariablespaths] = extended_path_core(ep.periods, M_.endo_nbr, M_.exo_nbr, innovations.positive_var_indx, ...
spfm_exo_simul, ep.init, endogenous_variables_paths(:,t-1), ...
oo_.steady_state, ...
verbosity, ep.stochastic.order, ...
M_, pfm, ep.stochastic.algo, ep.solve_algo, ep.stack_solve_algo, ...
options_.lmmcp, ...
options_, ...
oo_, initialguess);
if t>2
[endogenous_variables_paths(:,t), info_convergence, endogenousvariablespaths, y] = extended_path_core(innovations.positive_var_indx, ...
spfm_exo_simul, ...
endogenous_variables_paths(:,t-1), ...
pfm, ...
M_, ...
options_, ...
oo_, ...
initialguess, ...
y);
else
[endogenous_variables_paths(:,t), info_convergence, endogenousvariablespaths, y, pfm, options_] = extended_path_core(innovations.positive_var_indx, ...
spfm_exo_simul, ...
endogenous_variables_paths(:,t-1), ...
pfm, ...
M_, ...
options_, ...
oo_, ...
initialguess, ...
[]);
end
if ~info_convergence
msg = sprintf('No convergence of the (stochastic) perfect foresight solver (in period %s)!', int2str(t));
warning(msg)
......@@ -106,4 +123,4 @@ if any(isnan(endogenous_variables_paths(:)))
end
ts = dseries(transpose(endogenous_variables_paths), initial_period, M_.endo_names);
oo_.endo_simul = transpose(ts.data);
\ No newline at end of file
oo_.endo_simul = transpose(ts.data);
matlab/ep/extended_path_core.m
View file @ e9d4a40a
function [y, info_convergence, endogenousvariablespaths] = extended_path_core(periods,endo_nbr,exo_nbr,positive_var_indx, ...
exo_simul,init,initial_conditions,...
steady_state, ...
debug,order,M_,pfm,algo,solve_algo,stack_solve_algo,...
olmmcp,options_,oo_,initialguess)
function [y1, info_convergence, endogenousvariablespaths, y, pfm, options_] = extended_path_core(positive_var_indx, ...
exo_simul, ...
initial_conditions ,...
pfm, ...
M_, ...
options_, ...
oo_, ...
initialguess, ...
y)
% Copyright © 2016-2023 Dynare Team
% Copyright © 2016-2025 Dynare Team
%
% This file is part of Dynare.
%
......@@ -23,10 +27,22 @@ function [y, info_convergence, endogenousvariablespaths] = extended_path_core(pe
ep = options_.ep;
periods = ep.periods;
endo_nbr = M_.endo_nbr;
exo_nbr = M_.exo_nbr;
init = options_.ep.use_first_order_solution_as_initial_guess;
steady_state = oo_.steady_state;
debug = options_.verbosity;
order = options_.ep.stochastic.order;
algo = ep.stochastic.algo;
solve_algo = ep.solve_algo;
stack_solve_algo = ep.stack_solve_algo;
if init% Compute first order solution (Perturbation)...
endo_simul = simult_(M_,options_,initial_conditions,oo_.dr,exo_simul(2:end,:),1);
else
if nargin==19 && ~isempty(initialguess)
if nargin>7 && ~isempty(initialguess)
% Note that the first column of initialguess should be equal to initial_conditions.
endo_simul = initialguess;
else
......@@ -34,6 +50,10 @@ else
end
end
if nargin~=9
y = [];
end
oo_.endo_simul = endo_simul;
if debug
......@@ -41,30 +61,39 @@ if debug
end
if options_.bytecode && order > 0
error('Option order > 0 of extended_path command is not compatible with bytecode option.')
error('Option order > 0 of extended_path command with order>0 is not compatible with bytecode option.')
end
if options_.block && order > 0
error('Option order > 0 of extended_path command is not compatible with block option.')
error('Option order > 0 of extended_path command with order>0 is not compatible with block option.')
end
if order == 0
% Extended Path
options_.periods = periods;
options_.block = pfm.block;
oo_.endo_simul = endo_simul;
oo_.exo_simul = exo_simul;
oo_.steady_state = steady_state;
options_.lmmcp = olmmcp;
options_.solve_algo = solve_algo;
options_.stack_solve_algo = stack_solve_algo;
[endogenousvariablespaths, info_convergence] = perfect_foresight_solver_core(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
[endogenousvariablespaths, info_convergence] = perfect_foresight_solver_core(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, [], M_, options_);
else
% Stochastic Extended Path
switch(algo)
case 0
[flag, endogenousvariablespaths] = ...
solve_stochastic_perfect_foresight_model(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order);
case 1
[flag, endogenousvariablespaths] = ...
solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options_, pfm, ep.stochastic.order);
case 0
% Full tree of future trajectories.
if nargout>4
[flag, endogenousvariablespaths, errorcode, y, pfm, options_] = solve_stochastic_perfect_foresight_model_0(endo_simul, exo_simul, y, options_, M_, pfm);
else
[flag, endogenousvariablespaths, errorcode, y] = solve_stochastic_perfect_foresight_model_0(endo_simul, exo_simul, y, options_, M_, pfm);
end
case 1
% Sparse tree of future histories.
if nargout>4
[flag, endogenousvariablespaths, errorcode, y, pfm, options_] = solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, y, options_, M_, pfm);
else
[flag, endogenousvariablespaths, errorcode, y] = solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, y, options_, M_, pfm);
end
end
info_convergence = ~flag;
end
......@@ -74,7 +103,7 @@ if ~info_convergence && ~options_.no_homotopy
end
if info_convergence
y = endogenousvariablespaths(:,2);
y1 = endogenousvariablespaths(:,2);
else
y = NaN(size(endo_nbr,1));
y1 = NaN(size(endo_nbr,1));
end
matlab/ep/extended_path_homotopy.m
View file @ e9d4a40a
function [info_convergence, endo_simul] = extended_path_homotopy(endo_simul, exo_simul, M_, options_, oo_, pfm, ep, order, algo, method, debug)
% Copyright © 2016-2023 Dynare Team
% Copyright © 2016-2025 Dynare Team
%
% This file is part of Dynare.
%
......@@ -35,12 +35,12 @@ if ismember(method, [1, 2])
oo_.endo_simul(:,1) = oo_.steady_state + weight*(endo_simul0(:,1) - oo_.steady_state);
oo_.exo_simul = bsxfun(@plus, weight*exo_simul, (1-weight)*transpose(oo_.exo_steady_state));
if order==0
[endo_simul_new, success] = perfect_foresight_solver_core(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
[endo_simul_new, success] = perfect_foresight_solver_core(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, [], M_, options_);
else
switch(algo)
case 0
[flag, endo_simul_new] = ...
solve_stochastic_perfect_foresight_model(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order);
solve_stochastic_perfect_foresight_model_0(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order);
case 1
[flag, endo_simul_new] = ...
solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options_, pfm, ep.stochastic.order);
......@@ -102,12 +102,12 @@ if isequal(method, 3) || (isequal(method, 2) && noconvergence)
oo_.endo_simul = endo_simul;
oo_.exo_simul = bsxfun(@plus, weight*exo_simul, (1-weight)*transpose(oo_.exo_steady_state));
if order==0
[endo_simul_new, success] = perfect_foresight_solver_core(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
[endo_simul_new, success] = perfect_foresight_solver_core(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, [], M_, options_);
else
switch(algo)
case 0
[flag, endo_simul_new] = ...
solve_stochastic_perfect_foresight_model(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order);
solve_stochastic_perfect_foresight_model_0(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order);
case 1
[flag, endo_simul_new] = ...
solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options_, pfm, ep.stochastic.order);
......
matlab/ep/extended_path_initialization.m
View file @ e9d4a40a
function [initial_conditions, innovations, pfm, ep, verbosity, options_, oo_] = extended_path_initialization(initial_conditions, sample_size, exogenousvariables, options_, M_, oo_)
function [initial_conditions, innovations, pfm, options_, oo_] = extended_path_initialization(initial_conditions, sample_size, exogenousvariables, options_, M_, oo_)
% [initial_conditions, innovations, pfm, ep, verbosity, options_, oo_] = extended_path_initialization(initial_conditions, sample_size, exogenousvariables, options_, M_, oo_)
% Initialization of the extended path routines.
%
......@@ -16,7 +16,7 @@ function [initial_conditions, innovations, pfm, ep, verbosity, options_, oo_] =
%
% SPECIAL REQUIREMENTS
% Copyright © 2016-2024 Dynare Team
% Copyright © 2016-2025 Dynare Team
%
% This file is part of Dynare.
%
......@@ -36,8 +36,7 @@ function [initial_conditions, innovations, pfm, ep, verbosity, options_, oo_] =
ep = options_.ep;
% Set verbosity levels.
options_.verbosity = ep.verbosity;
verbosity = ep.verbosity+ep.debug;
options_.verbosity = ep.verbosity+ep.debug;
% Set maximum number of iterations for the deterministic solver.
options_.simul.maxit = ep.maxit;
......@@ -71,7 +70,7 @@ options_.stack_solve_algo = ep.stack_solve_algo;
% Compute the first order reduced form if needed.
dr = struct();
if ep.init
if ep.use_first_order_solution_as_initial_guess
options_.order = 1;
oo_.dr=set_state_space(dr,M_);
[oo_.dr,info,M_.params] = resol(0,M_,options_,oo_.dr,oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
......@@ -101,7 +100,11 @@ end
% hybrid correction
pfm.hybrid_order = ep.stochastic.hybrid_order;
if pfm.hybrid_order
if pfm.hybrid_order==1
warning('extended_path:: hybrid=1 is equivalent to hybrid=0 (option value must be an integer greater than 1 to be effective).')
end
if pfm.hybrid_order>1
oo_.dr = set_state_space(oo_.dr, M_);
options = options_;
options.order = pfm.hybrid_order;
......@@ -110,29 +113,37 @@ else
pfm.dr = [];
end
% Deactivate homotopy with SEP
if ep.stochastic.order>0
options_.no_homotopy = true;
end
% number of nonzero derivatives
pfm.nnzA = M_.NNZDerivatives(1);
% setting up integration nodes if order > 0
if ep.stochastic.order > 0
if ep.stochastic.order>0
[nodes,weights,nnodes] = setup_integration_nodes(options_.ep,pfm);
pfm.nodes = nodes;
pfm.weights = weights;
pfm.nnodes = nnodes;
% compute number of blocks
[block_nbr,pfm.world_nbr] = get_block_world_nbr(ep.stochastic.algo,nnodes,ep.stochastic.order,ep.periods);
[pfm.block_nbr, pfm.world_nbr] = get_block_world_nbr(ep.stochastic.algo,nnodes,ep.stochastic.order,ep.periods);
else
block_nbr = ep.periods;
end
% set boundaries if mcp
[lb, ub] = feval(sprintf('%s.dynamic_complementarity_conditions', M_.fname), M_.params);
pfm.eq_index = M_.dynamic_mcp_equations_reordering;
if options_.ep.solve_algo == 10
options_.lmmcp.lb = repmat(lb,block_nbr,1);
options_.lmmcp.ub = repmat(ub,block_nbr,1);
elseif options_.ep.solve_algo == 11
options_.mcppath.lb = repmat(lb,block_nbr,1);
options_.mcppath.ub = repmat(ub,block_nbr,1);
if ~ep.stochastic.order
[lb, ub] = feval(sprintf('%s.dynamic_complementarity_conditions', M_.fname), M_.params);
pfm.eq_index = M_.dynamic_mcp_equations_reordering;
if options_.ep.solve_algo == 10
options_.lmmcp.lb = repmat(lb, ep.periods, 1);
options_.lmmcp.ub = repmat(ub, ep.periods, 1);
elseif options_.ep.solve_algo == 11
options_.mcppath.lb = repmat(lb, ep.periods, 1);
options_.mcppath.ub = repmat(ub, ep.periods, 1);
end
else
% For SEP, boundaries are set in solve_stochastic_perfect_foresight_model_{0,1}.m
end
pfm.block_nbr = block_nbr;
matlab/ep/extended_path_mc.m
View file @ e9d4a40a
function Simulations = extended_path_mc(initialconditions, samplesize, replic, exogenousvariables, options_, M_, oo_)
% Simulations = extended_path_mc(initialconditions, samplesize, replic, exogenousvariables, options_, M_, oo_)
% Stochastic simulation of a non linear DSGE model using the Extended Path method (Fair and Taylor 1983). A time
% series of size T is obtained by solving T perfect foresight models.
%
......@@ -19,7 +19,7 @@ function Simulations = extended_path_mc(initialconditions, samplesize, replic, e
%
% SPECIAL REQUIREMENTS
% Copyright © 2016-2023 Dynare Team
% Copyright © 2016-2025 Dynare Team
%
% This file is part of Dynare.
%
......@@ -36,7 +36,7 @@ function Simulations = extended_path_mc(initialconditions, samplesize, replic, e
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
[initialconditions, innovations, pfm, ep, ~, options_, oo_] = ...
[initialconditions, innovations, pfm, options_, oo_] = ...
extended_path_initialization(initialconditions, samplesize, exogenousvariables, options_, M_, oo_);
% Check the dimension of the first input argument
......@@ -64,12 +64,12 @@ data = NaN(size(initialconditions, 1), samplesize+1, replic);
vexo = NaN(innovations.effective_number_of_shocks, samplesize+1, replic);
info = NaN(replic, 1);
if ep.parallel
if options_.ep.parallel
% Use the Parallel toolbox.
parfor i=1:replic
innovations_ = innovations;
oo__ = oo_;
[shocks, spfm_exo_simul, innovations_, oo__] = extended_path_shocks(innovations_, ep, exogenousvariables(:,:,i), samplesize, M_, options_, oo__);
[shocks, spfm_exo_simul, innovations_, oo__] = extended_path_shocks(innovations_, exogenousvariables(:,:,i), samplesize, M_, options_, oo__);
endogenous_variables_paths = NaN(M_.endo_nbr,samplesize+1);
endogenous_variables_paths(:,1) = initialconditions(:,1);
exogenous_variables_paths = NaN(innovations_.effective_number_of_shocks,samplesize+1);
......@@ -80,12 +80,21 @@ if ep.parallel
t = t+1;
spfm_exo_simul(2,:) = shocks(t-1,:);
exogenous_variables_paths(:,t) = shocks(t-1,:);
[endogenous_variables_paths(:,t), info_convergence] = extended_path_core(ep.periods, M_.endo_nbr, M_.exo_nbr, innovations_.positive_var_indx, ...
spfm_exo_simul, ep.init, endogenous_variables_paths(:,t-1), ...
oo__.steady_state, ...
ep.verbosity, ep.stochastic.order, ...
M_, pfm,ep.stochastic.algo, ep.solve_algo, ep.stack_solve_algo, ...
options_.lmmcp, options_, oo__);
if t>2
% Set initial guess for the solver (using the solution of the
% previous period problem).
initialguess = [endogenousvariablespaths(:, 2:end), oo_.steady_state];
else
initialguess = [];
end
[endogenous_variables_paths(:,t), info_convergence, endogenousvariablespaths] = extended_path_core(innovations.positive_var_indx, ...
spfm_exo_simul, ...
endogenous_variables_paths(:,t-1), ...
pfm, ...
M_, ...
options_, ...
oo__, ...
initialguess);
if ~info_convergence
msg = sprintf('No convergence of the (stochastic) perfect foresight solver (in period %s, iteration %s)!', int2str(t), int2str(i));
warning(msg)
......@@ -99,7 +108,7 @@ if ep.parallel
else
% Sequential approach.
for i=1:replic
[shocks, spfm_exo_simul, innovations, oo_] = extended_path_shocks(innovations, ep, exogenousvariables(:,:,i), samplesize, M_, options_, oo_);
[shocks, spfm_exo_simul, innovations, oo_] = extended_path_shocks(innovations, exogenousvariables(:,:,i), samplesize, M_, options_, oo_);
endogenous_variables_paths = NaN(M_.endo_nbr,samplesize+1);
endogenous_variables_paths(:,1) = initialconditions(:,1);
exogenous_variables_paths = NaN(innovations.effective_number_of_shocks,samplesize+1);
......@@ -109,12 +118,21 @@ else
t = t+1;
spfm_exo_simul(2,:) = shocks(t-1,:);
exogenous_variables_paths(:,t) = shocks(t-1,:);
[endogenous_variables_paths(:,t), info_convergence] = extended_path_core(ep.periods, M_.endo_nbr, M_.exo_nbr, innovations.positive_var_indx, ...
spfm_exo_simul, ep.init, endogenous_variables_paths(:,t-1), ...
oo_.steady_state, ...
ep.verbosity, ep.stochastic.order, ...
M_, pfm,ep.stochastic.algo, ep.solve_algo, ep.stack_solve_algo, ...
options_.lmmcp, options_, oo_);
if t>2
% Set initial guess for the solver (using the solution of the
% previous period problem).
initialguess = [endogenousvariablespaths(:, 2:end), oo_.steady_state];
else
initialguess = [];
end
[endogenous_variables_paths(:,t), info_convergence, endogenousvariablespaths] = extended_path_core(innovations.positive_var_indx, ...
spfm_exo_simul, ...
endogenous_variables_paths(:,t-1), ...
pfm, ...
M_, ...
options_, ...
oo_, ...
initialguess);
if ~info_convergence
msg = sprintf('No convergence of the (stochastic) perfect foresight solver (in period %s, iteration %s)!', int2str(t), int2str(i));
warning(msg)
......
matlab/ep/extended_path_shocks.m
View file @ e9d4a40a
function [shocks, spfm_exo_simul, innovations, oo_] = extended_path_shocks(innovations, ep, exogenousvariables, sample_size,M_,options_, oo_)
% [shocks, spfm_exo_simul, innovations, oo_] = extended_path_shocks(innovations, ep, exogenousvariables, sample_size,M_,options_, oo_)
% Copyright © 2016-2023 Dynare Team
function [shocks, spfm_exo_simul, innovations, oo_] = extended_path_shocks(innovations, exogenousvariables, sample_size, M_, options_, oo_)
% Copyright © 2016-2025 Dynare Team
%
% This file is part of Dynare.
%
......@@ -19,14 +19,14 @@ function [shocks, spfm_exo_simul, innovations, oo_] = extended_path_shocks(innov
% Simulate shocks.
if isempty(exogenousvariables)
switch ep.innovation_distribution
switch options_.ep.innovation_distribution
case 'gaussian'
shocks = zeros(sample_size, M_.exo_nbr);
shocks(:,innovations.positive_var_indx) = transpose(transpose(innovations.covariance_matrix_upper_cholesky)*randn(innovations.effective_number_of_shocks,sample_size));
case 'calibrated'
options = options_;
options.periods = options.ep.periods;
oo_local = make_ex_(M_,options,oo_);
oo_local = make_ex_(M_, options, oo_);
shocks = oo_local.exo_simul(2:end,:);
otherwise
error(['extended_path:: ' ep.innovation_distribution ' distribution for the structural innovations is not (yet) implemented!'])
......@@ -38,4 +38,4 @@ end
% Copy the shocks in exo_simul
oo_.exo_simul = shocks;
spfm_exo_simul = repmat(oo_.exo_steady_state',ep.periods+2,1);
\ No newline at end of file
spfm_exo_simul = repmat(oo_.exo_steady_state',options_.ep.periods+2,1);
matlab/ep/setup_stochastic_perfect_foresight_model_solver.m
View file @ e9d4a40a
......@@ -5,7 +5,7 @@ function pfm = setup_stochastic_perfect_foresight_model_solver(M_,options_,oo_)
% o options_ [struct] Dynare's options structure
% o oo_ [struct] Dynare's results structure
% Copyright © 2013-2023 Dynare Team
% Copyright © 2013-2025 Dynare Team
%
% This file is part of Dynare.
%
......
matlab/ep/solve_stochastic_perfect_foresight_model.m deleted 100644 → 0
View file @ 369cc0e9
function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model(endo_simul,exo_simul,pfm,nnodes,order)
% Copyright © 2012-2017 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/>.
flag = 0;
err = 0;
stop = 0;
params = pfm.params;
steady_state = pfm.steady_state;
ny = pfm.ny;
periods = pfm.periods;
dynamic_model = pfm.dynamic_model;
lead_lag_incidence = pfm.lead_lag_incidence;
lead_lag_incidence_t = transpose(lead_lag_incidence);
nyp = pfm.nyp;
nyf = pfm.nyf;
i_cols_1 = pfm.i_cols_1;
i_cols_j = pfm.i_cols_j;
i_cols_T = nonzeros(lead_lag_incidence(1:2,:)');
maxit = pfm.maxit_;
tolerance = pfm.tolerance;
verbose = pfm.verbose;
number_of_shocks = size(exo_simul,2);
[nodes,weights] = gauss_hermite_weights_and_nodes(nnodes);
if number_of_shocks>1
nodes = repmat(nodes,1,number_of_shocks)*chol(pfm.Sigma);
% to be fixed for Sigma ~= I
for i=number_of_shocks:-1:1
rr(i) = {nodes(:,i)};
ww(i) = {weights};
end
nodes = cartesian_product_of_sets(rr{:});
weights = prod(cartesian_product_of_sets(ww{:}),2);
nnodes = nnodes^number_of_shocks;
else
nodes = nodes*sqrt(pfm.Sigma);
end
if verbose
disp (' -----------------------------------------------------');
disp ('MODEL SIMULATION :');
fprintf('\n');
end
z = endo_simul(lead_lag_incidence_t(:)>0);
[~, jacobian] = dynamic_model(z, exo_simul, params,steady_state, 2);
% Each column of Y represents a different world
% The upper right cells are unused
% The first row block is ny x 1
% The second row block is ny x nnodes
% The third row block is ny x nnodes^2
% and so on until size ny x nnodes^order
world_nbr = nnodes^order;
Y = repmat(endo_simul(:),1,world_nbr);
% The columns of A map the elements of Y such that
% each block of Y with ny rows are unfolded column wise
dimension = ny*(sum(nnodes.^(0:order-1),2)+(periods-order)*world_nbr);
if order == 0
i_upd_r = (1:ny*periods);
i_upd_y = i_upd_r + ny;
else
i_upd_r = zeros(dimension,1);
i_upd_y = i_upd_r;
i_upd_r(1:ny) = (1:ny);
i_upd_y(1:ny) = ny+(1:ny);
i1 = ny+1;
i2 = 2*ny;
n1 = ny+1;
n2 = 2*ny;
for i=2:periods
for j=1:nnodes^min(i-1,order)
i_upd_r(i1:i2) = (n1:n2)+(j-1)*ny*periods;
i_upd_y(i1:i2) = (n1:n2)+ny+(j-1)*ny*(periods+2);
i1 = i2+1;
i2 = i2+ny;
end
n1 = n2+1;
n2 = n2+ny;
end
end
if rows(lead_lag_incidence)>2
icA = [find(lead_lag_incidence(1,:)) find(lead_lag_incidence(2,:))+world_nbr*ny ...
find(lead_lag_incidence(3,:))+2*world_nbr*ny]';
else
if nyf
icA = [find(lead_lag_incidence(2,:))+world_nbr*ny find(lead_lag_incidence(3,:))+2*world_nbr*ny ]';
else
icA = [find(lead_lag_incidence(1,:)) find(lead_lag_incidence(2,:))+world_nbr*ny ]';
end
end
h1 = clock;
for iter = 1:maxit
A1 = sparse([],[],[],ny*(sum(nnodes.^(0:order-1),2)+1),dimension,(order+1)*world_nbr*nnz(jacobian));
res = zeros(ny,periods,world_nbr);
i_rows = 1:ny;
i_cols = find(lead_lag_incidence');
i_cols_p = i_cols(1:nyp);
i_cols_s = i_cols(nyp+(1:ny));
i_cols_f = i_cols(nyp+ny+(1:nyf));
i_cols_Ap = i_cols_p;
i_cols_As = i_cols_s;
i_cols_Af = i_cols_f - ny;
for i = 1:order+1
i_w_p = 1;
for j = 1:nnodes^(i-1)
innovation = exo_simul;
if i > 1
innovation(i+1,:) = nodes(mod(j-1,nnodes)+1,:);
end
if i <= order
for k=1:nnodes
y = [Y(i_cols_p,i_w_p);
Y(i_cols_s,j);
Y(i_cols_f,(j-1)*nnodes+k)];
[d1,jacobian] = dynamic_model(y,innovation,params,steady_state,i+1);
if i == 1
% in first period we don't keep track of
% predetermined variables
i_cols_A = [i_cols_As - ny; i_cols_Af];
A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_1);
else
i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af];
A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_j);
end
res(:,i,j) = res(:,i,j)+weights(k)*d1;
i_cols_Af = i_cols_Af + ny;
end
else
y = [Y(i_cols_p,i_w_p);
Y(i_cols_s,j);
Y(i_cols_f,j)];
[d1,jacobian] = dynamic_model(y,innovation,params,steady_state,i+1);
if i == 1
% in first period we don't keep track of
% predetermined variables
i_cols_A = [i_cols_As - ny; i_cols_Af];
A1(i_rows,i_cols_A) = jacobian(:,i_cols_1);
else
i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af];
A1(i_rows,i_cols_A) = jacobian(:,i_cols_j);
end
res(:,i,j) = d1;
i_cols_Af = i_cols_Af + ny;
end
i_rows = i_rows + ny;
if mod(j,nnodes) == 0
i_w_p = i_w_p + 1;
end
if i > 1
if mod(j,nnodes) == 0
i_cols_Ap = i_cols_Ap + ny;
end
i_cols_As = i_cols_As + ny;
end
end
i_cols_p = i_cols_p + ny;
i_cols_s = i_cols_s + ny;
i_cols_f = i_cols_f + ny;
end
nzA = cell(periods,world_nbr);
for j=1:world_nbr
i_rows_y = find(lead_lag_incidence')+(order+1)*ny;
offset_c = ny*(sum(nnodes.^(0:order-1),2)+j-1);
offset_r = (j-1)*ny;
for i=order+2:periods
[d1,jacobian] = dynamic_model(Y(i_rows_y,j), ...
exo_simul,params, ...
steady_state,i+1);
if i == periods
[ir,ic,v] = find(jacobian(:,i_cols_T));
else
[ir,ic,v] = find(jacobian(:,i_cols_j));
end
nzA{i,j} = [offset_r+ir,offset_c+icA(ic), v]';
res(:,i,j) = d1;
i_rows_y = i_rows_y + ny;
offset_c = offset_c + world_nbr*ny;
offset_r = offset_r + world_nbr*ny;
end
end
err = max(abs(res(i_upd_r)));
if err < tolerance
stop = 1;
if verbose
fprintf('\n') ;
disp([' Total time of simulation :' num2str(etime(clock,h1))]) ;
fprintf('\n') ;
disp(' Convergency obtained.') ;
fprintf('\n') ;
end
flag = 0;% Convergency obtained.
endo_simul = reshape(Y(:,1),ny,periods+2);
break
end
A2 = [nzA{:}]';
A = [A1; sparse(A2(:,1),A2(:,2),A2(:,3),ny*(periods-order-1)*world_nbr,dimension)];
dy = -A\res(i_upd_r);
Y(i_upd_y) = Y(i_upd_y) + dy;
end
if ~stop
if verbose
fprintf('\n') ;
disp([' Total time of simulation :' num2str(etime(clock,h1))]) ;
fprintf('\n') ;
disp('WARNING : maximum number of iterations is reached (modify options_.simul.maxit).') ;
fprintf('\n') ;
end
flag = 1;% more iterations are needed.
endo_simul = 1;
end
if verbose
disp ('-----------------------------------------------------') ;
end
matlab/ep/solve_stochastic_perfect_foresight_model_0.m 0 → 100644
View file @ e9d4a40a
function [errorflag, endo_simul, errorcode, y, pfm, options_] = solve_stochastic_perfect_foresight_model_0(endo_simul, exo_simul, y, options_, M_, pfm)
% Copyright © 2025 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/>.
update_pfm_struct = false;
update_options_struct = false;
if nargout>4
update_pfm_struct = true;
end
if nargout>5
update_options_struct = true;
end
dynamic_model = pfm.dynamic_model;
periods = pfm.periods;
order = pfm.stochastic_order;
lead_lag_incidence = pfm.lead_lag_incidence;
lead_lag_incidence_t = transpose(lead_lag_incidence);
ny = pfm.ny;
nyp = pfm.nyp;
nyf = pfm.nyf;
i_cols_1 = pfm.i_cols_1;
i_cols_j = pfm.i_cols_j;
i_cols_T = nonzeros(lead_lag_incidence(1:2,:)');
nodes = pfm.nodes;
weights = pfm.weights;
nnodes = pfm.nnodes;
if update_pfm_struct
% make sure that there is a node equal to zero
% and permute nodes and weights to have zero first
k = find(sum(abs(nodes),2) < 1e-12);
if ~isempty(k)
nodes = [nodes(k,:); nodes(1:k-1,:); nodes(k+1:end,:)];
weights = [weights(k); weights(1:k-1); weights(k+1:end)];
else
error('there is no nodes equal to zero')
end
pfm.nodes = nodes;
pfm.weights = weights;
if pfm.hybrid_order>0
if pfm.hybrid_order==2
pfm.h_correction = 0.5*pfm.dr.ghs2(dr.inv_order_var);
elseif pfm.hybrid_order>2
pfm.h_correction = pfm.dr.g_0(pfm.dr.inv_order_var);
else
pfm.h_correction = 0;
end
else
pfm.h_correction = 0;
end
z = endo_simul(lead_lag_incidence_t(:)>0);
[~, jacobian] = dynamic_model(z, exo_simul, pfm.params, pfm.steady_state, 2);
world_nbr = nnodes^order;
% The columns of A map the elements of Y such that
% each block of Y with ny rows are unfolded column wise
dimension = ny*(sum(nnodes.^(0:order-1),2)+(periods-order)*world_nbr);
i_upd_r = zeros(dimension,1);
i_upd_y = i_upd_r;
i_upd_r(1:ny) = (1:ny);
i_upd_y(1:ny) = ny+(1:ny);
i1 = ny+1;
i2 = 2*ny;
n1 = ny+1;
n2 = 2*ny;
for i=2:periods
for j=1:nnodes^min(i-1,order)
i_upd_r(i1:i2) = (n1:n2)+(j-1)*ny*periods;
i_upd_y(i1:i2) = (n1:n2)+ny+(j-1)*ny*(periods+2);
i1 = i2+1;
i2 = i2+ny;
end
n1 = n2+1;
n2 = n2+ny;
end
if rows(lead_lag_incidence)>2
icA = [find(lead_lag_incidence(1,:)) find(lead_lag_incidence(2,:))+world_nbr*ny ...
find(lead_lag_incidence(3,:))+2*world_nbr*ny]';
else
if nyf
icA = [find(lead_lag_incidence(2,:))+world_nbr*ny find(lead_lag_incidence(3,:))+2*world_nbr*ny ]';
else
icA = [find(lead_lag_incidence(1,:)) find(lead_lag_incidence(2,:))+world_nbr*ny ]';
end
end
i_rows = 1:ny;
i_cols = find(lead_lag_incidence');
i_cols_p = i_cols(1:nyp);
i_cols_s = i_cols(nyp+(1:ny));
i_cols_f = i_cols(nyp+ny+(1:nyf));
pfm.i_cols_Ap = i_cols_p;
pfm.i_cols_As = i_cols_s;
pfm.i_cols_Af = i_cols_f - ny;
pfm.i_hc = i_cols_f - 2*ny;
pfm.i_cols_p = i_cols_p;
pfm.i_cols_s = i_cols_s;
pfm.i_cols_f = i_cols_f;
pfm.i_rows = i_rows;
pfm.A1 = sparse([],[],[],ny*(sum(nnodes.^(0:order-1),2)+1),dimension,(order+1)*world_nbr*nnz(jacobian));
pfm.res = zeros(ny,periods,world_nbr);
pfm.order = order;
pfm.world_nbr = world_nbr;
pfm.i_cols_1 = i_cols_1;
pfm.i_cols_h = i_cols_j;
pfm.icA = icA;
pfm.i_cols_T = i_cols_T;
pfm.i_upd_r = i_upd_r;
pfm.i_upd_y = i_upd_y;
pfm.dimension = dimension;
end
pfm.Y = repmat(endo_simul(:),1,pfm.world_nbr);
if isempty(y)
y = repmat(pfm.steady_state, pfm.dimension/pfm.ny, 1);
end
if update_options_struct
% Set algorithm
options_.solve_algo = options_.ep.solve_algo;
options_.simul.maxit = options_.ep.maxit;
[lb, ub] = feval(sprintf('%s.dynamic_complementarity_conditions', M_.fname), pfm.params);
pfm.eq_index = M_.dynamic_mcp_equations_reordering;
if options_.ep.solve_algo == 10
options_.lmmcp.lb = repmat(lb, pfm.dimension/pfm.ny, 1);
options_.lmmcp.ub = repmat(ub, pfm.dimension/pfm.ny, 1);
elseif options_.ep.solve_algo == 11
options_.mcppath.lb = repmat(lb, pfm.dimension/pfm.ny, 1);
options_.mcppath.ub = repmat(ub, pfm.dimension/pfm.ny, 1);
end
end
[y, errorflag, ~, ~, errorcode] = dynare_solve(@ep_problem_0, y, options_.simul.maxit, options_.dynatol.f, options_.dynatol.x, options_, exo_simul, pfm);
endo_simul(:,2) = y(1:ny);
matlab/ep/solve_stochastic_perfect_foresight_model_1.m
View file @ e9d4a40a
function [errorflag, endo_simul, errorcode, y] = solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, Options, pfm, order, varargin)
function [errorflag, endo_simul, errorcode, y, pfm, options_] = solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, y, options_, M_, pfm)
% Copyright © 2012-2022 Dynare Team
% Copyright © 2012-2025 Dynare Team
%
% This file is part of Dynare.
%
......@@ -17,87 +17,74 @@ function [errorflag, endo_simul, errorcode, y] = solve_stochastic_perfect_foresi
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
if nargin < 6
homotopy_parameter = 1;
else
homotopy_parameter = varargin{1};
update_pfm_struct = false;
update_options_struct = false;
if nargout>4
update_pfm_struct = true;
end
flag = 0;
err = 0;
EpOptions = Options.ep;
params = pfm.params;
steady_state = pfm.steady_state;
ny = pfm.ny;
periods = pfm.periods;
dynamic_model = pfm.dynamic_model;
lead_lag_incidence = pfm.lead_lag_incidence;
nyp = pfm.nyp;
nyf = pfm.nyf;
i_cols_1 = pfm.i_cols_1;
i_cols_A1 = pfm.i_cols_A1;
i_cols_j = pfm.i_cols_j;
i_cols_T = nonzeros(lead_lag_incidence(1:2,:)');
hybrid_order = pfm.hybrid_order;
dr = pfm.dr;
nodes = pfm.nodes;
weights = pfm.weights;
nnodes = pfm.nnodes;
maxit = pfm.maxit_;
tolerance = pfm.tolerance;
verbose = pfm.verbose;
number_of_shocks = size(exo_simul,2);
% make sure that there is a node equal to zero
% and permute nodes and weights to have zero first
k = find(sum(abs(nodes),2) < 1e-12);
if ~isempty(k)
nodes = [nodes(k,:); nodes(1:k-1,:); nodes(k+1:end,:)];
weights = [weights(k); weights(1:k-1); weights(k+1:end)];
else
error('there is no nodes equal to zero')
if nargout>5
update_options_struct = true;
end
if hybrid_order > 0
if hybrid_order == 2
h_correction = 0.5*dr.ghs2(dr.inv_order_var);
if update_pfm_struct
periods = pfm.periods;
order = pfm.stochastic_order;
ny = pfm.ny;
lead_lag_incidence = pfm.lead_lag_incidence;
i_cols_1 = pfm.i_cols_1;
i_cols_j = pfm.i_cols_j;
i_cols_T = nonzeros(lead_lag_incidence(1:2,:)');
nodes = pfm.nodes;
weights = pfm.weights;
nnodes = pfm.nnodes;
% Make sure that there is a node equal to zero and permute nodes and weights to have zero first
k = find(sum(abs(nodes),2) < 1e-12);
if ~isempty(k)
nodes = [nodes(k,:); nodes(1:k-1,:); nodes(k+1:end,:)];
weights = [weights(k); weights(1:k-1); weights(k+1:end)];
else
error('there is no nodes equal to zero')
end
else
h_correction = 0;
end
if verbose
disp (' -----------------------------------------------------');
disp ('MODEL SIMULATION :');
fprintf('\n');
end
pfm.nodes = nodes;
pfm.weights = weights;
% Each column of Y represents a different world
% The upper right cells are unused
% The first row block is ny x 1
% The second row block is ny x nnodes
% The third row block is ny x nnodes^2
% and so on until size ny x nnodes^order
world_nbr = pfm.world_nbr;
Y = endo_simul(:,2:end-1);
Y = repmat(Y,1,world_nbr);
pfm.y0 = endo_simul(:,1);
if pfm.hybrid_order > 0
if pfm.hybrid_order == 2
pfm.h_correction = 0.5*pfm.dr.ghs2(pfm.dr.inv_order_var);
elseif pfm.hybrid_order>2
pfm.h_correction = pfm.dr.g_0(pfm.dr.inv_order_var);
else
pfm.h_correction = 0;
end
else
pfm.h_correction = 0;
end
% Each column of Y represents a different world
% The upper right cells are unused
% The first row block is ny x 1
% The second row block is ny x nnodes
% The third row block is ny x nnodes^2
% and so on until size ny x nnodes^order
world_nbr = pfm.world_nbr;
% The columns of A map the elements of Y such that
% each block of Y with ny rows are unfolded column wise
% number of blocks
block_nbr = pfm.block_nbr;
% dimension of the problem
dimension = ny*block_nbr;
pfm.dimension = dimension;
if order == 0
i_upd_r = (1:ny*periods)';
i_upd_y = i_upd_r + ny;
else
i_upd_r = zeros(dimension,1);
% The columns of A map the elements of Y such that
% each block of Y with ny rows are unfolded column wise
% number of blocks
block_nbr = pfm.block_nbr;
% dimension of the problem
dimension = ny*block_nbr;
pfm.dimension = dimension;
i_upd_r = zeros(dimension, 1);
i_upd_y = i_upd_r;
i_upd_r(1:ny) = (1:ny);
i_upd_y(1:ny) = ny+(1:ny);
......@@ -116,32 +103,40 @@ else
n1 = n2+1;
n2 = n2+ny;
end
icA = [find(lead_lag_incidence(1,:)) find(lead_lag_incidence(2,:))+world_nbr*ny ...
find(lead_lag_incidence(3,:))+2*world_nbr*ny]';
pfm.i_rows = 1:ny;
pfm.i_cols = find(lead_lag_incidence');
pfm.i_cols_1 = i_cols_1;
pfm.i_cols_h = i_cols_j;
pfm.icA = icA;
pfm.i_cols_T = i_cols_T;
pfm.i_upd_r = i_upd_r;
pfm.i_upd_y = i_upd_y;
end
if isempty(y)
y = repmat(pfm.steady_state, pfm.block_nbr, 1);
end
if update_options_struct
% Set algorithm
options_.solve_algo = options_.ep.solve_algo;
options_.simul.maxit = options_.ep.maxit;
[lb, ub] = feval(sprintf('%s.dynamic_complementarity_conditions', M_.fname), pfm.params);
pfm.eq_index = M_.dynamic_mcp_equations_reordering;
if options_.ep.solve_algo == 10
options_.lmmcp.lb = repmat(lb, pfm.block_nbr, 1);
options_.lmmcp.ub = repmat(ub, pfm.block_nbr, 1);
elseif options_.ep.solve_algo == 11
options_.mcppath.lb = repmat(lb, pfm.block_nbr, 1);
options_.mcppath.ub = repmat(ub, pfm.block_nbr, 1);
end
end
icA = [find(lead_lag_incidence(1,:)) find(lead_lag_incidence(2,:))+world_nbr*ny ...
find(lead_lag_incidence(3,:))+2*world_nbr*ny]';
h1 = clock;
pfm.order = order;
pfm.world_nbr = world_nbr;
pfm.nodes = nodes;
pfm.nnodes = nnodes;
pfm.weights = weights;
pfm.h_correction = h_correction;
pfm.i_rows = 1:ny;
i_cols = find(lead_lag_incidence');
pfm.i_cols = i_cols;
pfm.nyp = nyp;
pfm.nyf = nyf;
pfm.hybrid_order = hybrid_order;
pfm.i_cols_1 = i_cols_1;
pfm.i_cols_h = i_cols_j;
pfm.icA = icA;
pfm.i_cols_T = i_cols_T;
pfm.i_upd_r = i_upd_r;
pfm.i_upd_y = i_upd_y;
Options.steady.maxit = 100;
y = repmat(steady_state,block_nbr,1);
Options.solve_algo = Options.ep.solve_algo;
Options.steady.maxit = Options.ep.maxit;
[y, errorflag, ~, ~, errorcode] = dynare_solve(@ep_problem_2, y, Options.simul.maxit, Options.dynatol.f, Options.dynatol.x, Options, exo_simul, pfm);
endo_simul(:,2) = y(1:ny);
\ No newline at end of file
pfm.y0 = endo_simul(:,1);
[y, errorflag, ~, ~, errorcode] = dynare_solve(@ep_problem_1, y, options_.simul.maxit, options_.dynatol.f, options_.dynatol.x, options_, exo_simul, pfm);
endo_simul(:,2) = y(1:pfm.ny);
matlab/estimation/GetAllPosteriorDraws.m
View file @ e9d4a40a
......@@ -51,6 +51,18 @@ if issmc(options_)
else
draws = posterior.tlogpostkernel;
end
elseif isonline(options_)
% Load draws from the posterior distribution
posterior = load(sprintf('%s/online/parameters_particles_final.mat', dname));
if column>0 || strcmp(column,'all')
if strcmp(column,'all')
draws = transpose(posterior.param);
else
draws = transpose(posterior.param(column,:));
end
else
draws=NaN(size(posterior.param,2),1);
end
elseif isdime(options_)
posterior = load(sprintf('%s%s%s%schains.mat', dname, filesep(), 'dime', filesep()));
tune = posterior.tune;
......
matlab/estimation/GetPosteriorMeanVariance.m
View file @ e9d4a40a
......@@ -36,6 +36,12 @@ if issmc(options_)
mean = sum(posterior.particles, 2)/length(posterior.tlogpostkernel);
% Compute the posterior covariance
variance = (posterior.particles-mean)*(posterior.particles-mean)'/length(posterior.tlogpostkernel);
elseif isonline(options_)
posterior = load(sprintf('%s/online/parameters_particles_final.mat', M_.dname));
% Compute the posterior mean
mean = sum(posterior.param, 2)/size(posterior.param, 2);
% Compute the posterior covariance
variance = (posterior.param-mean)*(posterior.param-mean)'/size(posterior.param, 2);
elseif isdime(options_)
posterior = load(sprintf('%s%s%s%schains.mat', M_.dname, filesep(), 'dime', filesep()));
tune = posterior.tune;
......