diff --git a/matlab/ep/ep_problem_1.m b/matlab/ep/ep_problem_1.m
new file mode 100644
index 0000000000000000000000000000000000000000..51c7fb1a80e8b955b90629151913c26a082b9fa8
--- /dev/null
+++ b/matlab/ep/ep_problem_1.m
@@ -0,0 +1,122 @@
+function [res, A, info] = ep_problem_1(y, x, pfm)
+
+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==2 && 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);
diff --git a/matlab/ep/ep_problem_2.m b/matlab/ep/ep_problem_2.m
index 9168631bfe5a06927da76bb836f360d2357e36b6..d94f277f8b374810675d9c8bfda99fafb43ac342 100644
--- a/matlab/ep/ep_problem_2.m
+++ b/matlab/ep/ep_problem_2.m
@@ -1,37 +1,37 @@
-function [res,A,info] = ep_problem_2(y,x,pm)
+function [res,A,info] = ep_problem_2(y, x, pfm)
-info = 0;
-res = [];
+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;
@@ -39,10 +39,11 @@ 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 = pfm.Y; %zeros(ny*(periods+2),world_nbr);
+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;
@@ -81,7 +82,7 @@ for i = 1:order+1
if hybrid_order == 2 && (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 +193,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);
diff --git a/matlab/ep/extended_path.m b/matlab/ep/extended_path.m
index 2b807d36f8b3de236ab821cbd7284f85cac94f1e..c4a590d0cb749aa5da5c5e458293c271dd9f1653 100644
--- a/matlab/ep/extended_path.m
+++ b/matlab/ep/extended_path.m
@@ -36,18 +36,18 @@ 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_] = ...
+[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 +67,25 @@ 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.use_first_order_solution_as_initial_guess, 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] = extended_path_core(innovations.positive_var_indx, ...
+ spfm_exo_simul, ...
+ endogenous_variables_paths(:,t-1), ...
+ pfm, ...
+ M_, ...
+ options_, ...
+ oo_, ...
+ initialguess);
+ else
+ [endogenous_variables_paths(:,t), info_convergence, endogenousvariablespaths, 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)
diff --git a/matlab/ep/extended_path_core.m b/matlab/ep/extended_path_core.m
index f5bf14d491250de9964d9a7a8258d5a3e0266abf..71571b112c14f9835574f2e6f6ad8e770ce40c55 100644
--- a/matlab/ep/extended_path_core.m
+++ b/matlab/ep/extended_path_core.m
@@ -1,10 +1,13 @@
-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 [y, info_convergence, endogenousvariablespaths, pfm, options_] = extended_path_core(positive_var_indx, ...
+ exo_simul, ...
+ initial_conditions ,...
+ pfm, ...
+ M_, ...
+ options_, ...
+ oo_, ...
+ initialguess)
-% Copyright © 2016-2023 Dynare Team
+% Copyright © 2016-2025 Dynare Team
%
% This file is part of Dynare.
%
@@ -23,6 +26,18 @@ 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
@@ -41,30 +56,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_);
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>3
+ [flag, endogenousvariablespaths, errorcode, ~, pfm, options_] = solve_stochastic_perfect_foresight_model_0(endo_simul, exo_simul, options_, M_, pfm);
+ else
+ [flag, endogenousvariablespaths, errorcode] = solve_stochastic_perfect_foresight_model_0(endo_simul, exo_simul, options_, M_, pfm);
+ end
+ case 1
+ % Sparse tree of future histories.
+ if nargout>3
+ [flag, endogenousvariablespaths, errorcode, ~, pfm, options_] = solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options_, M_, pfm);
+ else
+ [flag, endogenousvariablespaths, errorcode] = solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options_, M_, pfm);
+ end
end
info_convergence = ~flag;
end
diff --git a/matlab/ep/extended_path_initialization.m b/matlab/ep/extended_path_initialization.m
index 0359b9ee23ad5f1f61e4e29ebb5271def398ba89..90828dc4b6895472bf0af0170f4dae43b2e37d05 100644
--- a/matlab/ep/extended_path_initialization.m
+++ b/matlab/ep/extended_path_initialization.m
@@ -1,4 +1,4 @@
-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.
%
@@ -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;
@@ -114,25 +113,28 @@ end
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;
diff --git a/matlab/ep/extended_path_mc.m b/matlab/ep/extended_path_mc.m
index b49f4bcbba5c3ca689a5d392352b08a5204c0845..8fbf172563e6602ad31960060f073621c1b1c2d1 100644
--- a/matlab/ep/extended_path_mc.m
+++ b/matlab/ep/extended_path_mc.m
@@ -1,5 +1,5 @@
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.
%
@@ -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.use_first_order_solution_as_initial_guess, 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.use_first_order_solution_as_initial_guess, 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)
diff --git a/matlab/ep/extended_path_shocks.m b/matlab/ep/extended_path_shocks.m
index 7c567ee4e1e3f20bf474277b71aa50153e5748bc..f86d89c300b452ab4b7d9b573797ef3e3dfb7459 100644
--- a/matlab/ep/extended_path_shocks.m
+++ b/matlab/ep/extended_path_shocks.m
@@ -1,6 +1,6 @@
-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);
diff --git a/matlab/ep/setup_stochastic_perfect_foresight_model_solver.m b/matlab/ep/setup_stochastic_perfect_foresight_model_solver.m
index 5266a1b4f0ab9b457a590ff962505d5af9eca255..b319672274e5f77e99b7fb2813270cda9f950882 100644
--- a/matlab/ep/setup_stochastic_perfect_foresight_model_solver.m
+++ b/matlab/ep/setup_stochastic_perfect_foresight_model_solver.m
@@ -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.
%
diff --git a/matlab/ep/solve_stochastic_perfect_foresight_model.m b/matlab/ep/solve_stochastic_perfect_foresight_model.m
index d18b6cf05b94a32fed1eacb375e24a481b852785..8263f5e706c3a2f6c49a264817dba545858cf722 100644
--- a/matlab/ep/solve_stochastic_perfect_foresight_model.m
+++ b/matlab/ep/solve_stochastic_perfect_foresight_model.m
@@ -33,6 +33,16 @@ 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,:)');
+hybrid_order = pfm.hybrid_order;
+if hybrid_order > 0
+ if hybrid_order == 2
+ h_correction = 0.5*pfm.dr.ghs2(pfm.dr.inv_order_var);
+ end
+else
+ h_correction = 0;
+end
+%h_correction = pfm.h_correction;
+
maxit = pfm.maxit_;
tolerance = pfm.tolerance;
@@ -122,6 +132,7 @@ for iter = 1:maxit
i_cols_Ap = i_cols_p;
i_cols_As = i_cols_s;
i_cols_Af = i_cols_f - ny;
+ i_hc = i_cols_f - 2*ny;
for i = 1:order+1
i_w_p = 1;
for j = 1:nnodes^(i-1)
@@ -148,9 +159,16 @@ for iter = 1:maxit
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)];
+ % i==order+1
+ if hybrid_order==2
+ y = [Y(i_cols_p,i_w_p);
+ Y(i_cols_s,j);
+ Y(i_cols_f,j)+h_correction(i_hc)];
+ else
+ y = [Y(i_cols_p,i_w_p);
+ Y(i_cols_s,j);
+ Y(i_cols_f,j)];
+ end
[d1,jacobian] = dynamic_model(y,innovation,params,steady_state,i+1);
if i == 1
% in first period we don't keep track of
diff --git a/matlab/ep/solve_stochastic_perfect_foresight_model_0.m b/matlab/ep/solve_stochastic_perfect_foresight_model_0.m
new file mode 100644
index 0000000000000000000000000000000000000000..02947bb21b63db46055a9b2f739529c9ef50f990
--- /dev/null
+++ b/matlab/ep/solve_stochastic_perfect_foresight_model_0.m
@@ -0,0 +1,170 @@
+function [errorflag, endo_simul, errorcode, y, pfm, options_] = solve_stochastic_perfect_foresight_model_0(endo_simul, exo_simul, 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;
+
+ hybrid_order = pfm.hybrid_order;
+ dr = pfm.dr;
+ if hybrid_order > 0
+ if hybrid_order == 2
+ h_correction = 0.5*dr.ghs2(dr.inv_order_var);
+ end
+ else
+ h_correction = 0;
+ end
+
+ pfm.h_correction = h_correction;
+
+ z = endo_simul(lead_lag_incidence_t(:)>0);
+ [~, jacobian] = dynamic_model(z, exo_simul, pfm.params, pfm.steady_state, 2);
+
+ 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);
+
+ 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.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;
+ pfm.Y = Y;
+
+ pfm.dimension = dimension;
+
+end
+
+y = repmat(pfm.steady_state, pfm.dimension/pfm.ny, 1);
+
+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_1, y, options_.simul.maxit, options_.dynatol.f, options_.dynatol.x, options_, exo_simul, pfm);
+
+endo_simul(:,2) = y(1:ny);
diff --git a/matlab/ep/solve_stochastic_perfect_foresight_model_1.m b/matlab/ep/solve_stochastic_perfect_foresight_model_1.m
index 633cd542c6347a60449211fea13f5d8fd0379aff..d65c464757f618f76f8c8e097334869d8620068e 100644
--- a/matlab/ep/solve_stochastic_perfect_foresight_model_1.m
+++ b/matlab/ep/solve_stochastic_perfect_foresight_model_1.m
@@ -1,6 +1,6 @@
-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, options_, M_, pfm)
-% Copyright © 2012-2022 Dynare Team
+% Copyright © 2012-2025 Dynare Team
%
% This file is part of Dynare.
%
@@ -17,87 +17,73 @@ 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);
-
-% 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);
+ if pfm.hybrid_order > 0
+ if pfm.hybrid_order == 2
+ pfm.h_correction = 0.5*pfm.dr.ghs2(pfm.dr.inv_order_var);
+ 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;
+ % Y = endo_simul(:,2:end-1);
+ pfm.Y = repmat(endo_simul(:), 1, world_nbr);
+ pfm.y0 = endo_simul(:,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 +102,36 @@ 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
+
+y = repmat(pfm.steady_state, pfm.block_nbr, 1);
+
+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
+
+[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:pfm.ny);
diff --git a/tests/ep/burnside.mod b/tests/ep/burnside.mod
index 958cb8f62f1c17ddc9087ba84dc3d22dd5ac92cb..a932a95335bb500b36a3db138b175bfc76307f5e 100644
--- a/tests/ep/burnside.mod
+++ b/tests/ep/burnside.mod
@@ -34,6 +34,7 @@ end
set_dynare_seed('default');
stoch_simul(order=1,irf=0,periods=201);
y_perturbation_1 = oo_.endo_simul(1,:)';
+e_1 = oo_.exo_simul;
set_dynare_seed('default');
stoch_simul(order=2,irf=0,periods=201);
@@ -48,16 +49,31 @@ options_.simul.maxit = 100;
set_dynare_seed('default');
ytrue=exact_solution(M_,oo_, 800);
-
+/*
set_dynare_seed('default');
options_.ep.stochastic.order = 0;
-ts = extended_path([], 200, [], options_, M_, oo_);
+ts = extended_path([], 200, e_1, options_, M_, oo_);
+*/
+
set_dynare_seed('default');
-options_.ep.stochastic.order = 2;
-options_.ep.stochastic.IntegrationAlgorithm='Stroud-Cubature-3';//'Unscented'; //'Tensor-Gaussian-Quadrature';
+options_.ep.stochastic.order = 1;
+options_.ep.stochastic.IntegrationAlgorithm='Tensor-Gaussian-Quadrature';//'Stroud-Cubature-3';//'Unscented'; //;
+options_.ep.stochastic.quadrature.nodes = 3;
+options_.ep.stochastic.algo = 0;
+options_.ep.stochastic.hybrid_order = 2;
+ts1 = extended_path([], 200, e_1, options_, M_, oo_);
+
+
+set_dynare_seed('default');
+options_.ep.stochastic.order = 1;
+options_.ep.stochastic.IntegrationAlgorithm='Tensor-Gaussian-Quadrature';//'Stroud-Cubature-3';//'Unscented'; //'Tensor-Gaussian-Quadrature';
options_.ep.stochastic.quadrature.nodes = 3;
-ts1_4 = extended_path([], 200, [], options_, M_, oo_);
+options_.ep.stochastic.algo = 1;
+options_.ep.stochastic.hybrid_order = 2;
+ts1_ = extended_path([], 200, e_1, options_, M_, oo_);
+
+return
set_dynare_seed('default');
options_.ep.stochastic.order = 4;
diff --git a/tests/ep/rbc.mod b/tests/ep/rbc.mod
index 4003c4aa0281499b6b8dfc9e3990ddf4c94597e7..c1b29b77b28b91043de5ca09656fab74f3923c1d 100644
--- a/tests/ep/rbc.mod
+++ b/tests/ep/rbc.mod
@@ -72,12 +72,13 @@ end;
steady(nocheck);
-options_.ep.verbosity = 0;
options_.ep.stochastic.order = 0;
+options_.ep.stack_solve_algo=0;
ts0 = extended_path([], 10, [], options_, M_, oo_);
options_.ep.stochastic.order = 1;
+
options_.ep.stochastic.nodes = 3;
options_.ep.stochastic.IntegrationAlgorithm='Tensor-Gaussian-Quadrature';
ts1_3 = extended_path([], 10, [], options_, M_, oo_);
diff --git a/tests/ep/rbc_mc.mod b/tests/ep/rbc_mc.mod
index 5f30d6ce5e92d64016472bdd5e712929e3de3eaf..a0dcadbb2b15c42c943c49a6554474b798c82284 100644
--- a/tests/ep/rbc_mc.mod
+++ b/tests/ep/rbc_mc.mod
@@ -72,4 +72,6 @@ end;
steady(nocheck);
-Simulations = extended_path_mc([], 2, 2, [], options_, M_, oo_);
\ No newline at end of file
+options_.ep.stack_solve_algo=0;
+
+Simulations = extended_path_mc([], 10, 2, [], options_, M_, oo_);
diff --git a/tests/ep/rbcii.mod b/tests/ep/rbcii.mod
index 3f3fa32a5ad1d46c75281b34ffea87020efd8bb3..24385ccf196505b51f150d5b78f0f42b3fc9d8d0 100644
--- a/tests/ep/rbcii.mod
+++ b/tests/ep/rbcii.mod
@@ -50,13 +50,10 @@ end;
steady(nocheck);
- options_.ep.stochastic.order = 0;
-
- ts = extended_path([], 200, [], options_, M_, oo_);
- ts.save('rbcii-sim-data');
-
- options_.ep.stochastic.order = 1;
- ts1_4 = extended_path([], 200, [], options_, M_, oo_);
+ //options_.ep.stochastic.order = 0;
+ set_dynare_seed(2009);
+ Simulated_time_series = extended_path([], 200, [], options_, M_, oo_);
+ Simulated_time_series.save('rbcii-sim-data');
@#else
diff --git a/tests/ep/rbcii_MCP.mod b/tests/ep/rbcii_MCP.mod
index ffad81e92673981cc5e54ab7730ae8b676a8c958..01ec2854075673d2f1fe10aec9da64bb7ccdfea7 100644
--- a/tests/ep/rbcii_MCP.mod
+++ b/tests/ep/rbcii_MCP.mod
@@ -49,9 +49,10 @@ shocks;
stderr 0.10;
end;
-extended_path(periods=200,lmmcp);
+set_dynare_seed(2009);
+extended_path(periods=200,order=0,lmmcp);
-if any(oo_.endo_simul(strmatch('i',M_.endo_names,'exact'),:)<-1e-6)
+if any(Simulated_time_series.i.data<-1e-6)
error('lmmcp tag did not work.')
end
@@ -59,10 +60,10 @@ ds = dseries('rbcii-sim-data.mat');
if isoctave
tolerance=5e-5;
else
- tolerance=1e-6;
+ tolerance=1e-5;
end
-if any(abs(transpose(oo_.endo_simul(strmatch('i',M_.endo_names,'exact'),:))-ds.Investment.data)>tolerance)
+if any(abs(Simulated_time_series.i.data-ds.Investment.data)>tolerance)
error('Simulation with lmmcp returns different results.')
end