diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
index 8e38db33f4352d939c5c90eca518527c4890e67b..774754bd26a0a0b897953b0b88a35d5c504e1cf5 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
@@ -86,72 +86,12 @@ else
[oo_.endo_simul, oo_.deterministic_simulation] = ...
sim1_lbj(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
elseif options_.stack_solve_algo == 7
- periods = options_.periods;
- if ~isfield(options_.lmmcp,'lb')
- [lb,ub,pfm.eq_index] = get_complementarity_conditions(M_,options_.ramsey_policy);
- options_.lmmcp.lb = repmat(lb,periods,1);
- options_.lmmcp.ub = repmat(ub,periods,1);
- end
- y = oo_.endo_simul;
- y0 = y(:,1);
- yT = y(:,periods+2);
- z = y(:,2:periods+1);
- illi = M_.lead_lag_incidence';
- [i_cols,junk,i_cols_j] = find(illi(:));
- illi = illi(:,2:3);
- [i_cols_J1,junk,i_cols_1] = find(illi(:));
- i_cols_T = nonzeros(M_.lead_lag_incidence(1:2,:)');
- if options_.linear_approximation
- y_steady_state = oo_.steady_state;
- x_steady_state = transpose(oo_.exo_steady_state);
- ip = find(M_.lead_lag_incidence(1,:)');
- ic = find(M_.lead_lag_incidence(2,:)');
- in = find(M_.lead_lag_incidence(3,:)');
- % Evaluate the Jacobian of the dynamic model at the deterministic steady state.
- model_dynamic = str2func([M_.fname,'_dynamic']);
- [d1,jacobian] = model_dynamic(y_steady_state([ip; ic; in]), x_steady_state, M_.params, y_steady_state, 1);
- % Check that the dynamic model was evaluated at the steady state.
- if max(abs(d1))>1e-12
- error('Jacobian is not evaluated at the steady state!')
- end
- nyp = nnz(M_.lead_lag_incidence(1,:)) ;
- ny0 = nnz(M_.lead_lag_incidence(2,:)) ;
- nyf = nnz(M_.lead_lag_incidence(3,:)) ;
- nd = nyp+ny0+nyf; % size of y (first argument passed to the dynamic file).
- jexog = transpose(nd+(1:M_.exo_nbr));
- jendo = transpose(1:nd);
- z = bsxfun(@minus,z,y_steady_state);
- x = bsxfun(@minus,oo_.exo_simul,x_steady_state);
- [y,info] = dynare_solve(@linear_perfect_foresight_problem,z(:), options_, ...
- jacobian, y0-y_steady_state, yT-y_steady_state, ...
- x, M_.params, y_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_.NNZDerivatives(1),jendo,jexog);
- else
- [y,info] = dynare_solve(@perfect_foresight_problem,z(:),options_, ...
- 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_.NNZDerivatives(1));
- end
- if all(imag(y)<.1*options_.dynatol.f)
- if ~isreal(y)
- y = real(y);
- end
- else
- info = 1;
- end
if options_.linear_approximation
- oo_.endo_simul = [y0 bsxfun(@plus,reshape(y,M_.endo_nbr,periods),y_steady_state) yT];
- else
- oo_.endo_simul = [y0 reshape(y,M_.endo_nbr,periods) yT];
- end
- if info == 1
- oo_.deterministic_simulation.status = false;
+ [oo_.endo_simul, oo_.deterministic_simulation] = ...
+ solve_stacked_linear_problem(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
else
- oo_.deterministic_simulation.status = true;
+ [oo_.endo_simul, oo_.deterministic_simulation] = ...
+ solve_stacked_problem(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
end
end
end
diff --git a/matlab/perfect-foresight-models/private/initialize_stacked_problem.m b/matlab/perfect-foresight-models/private/initialize_stacked_problem.m
new file mode 100644
index 0000000000000000000000000000000000000000..d4ab762beba0108b527fbbcabe3174e04ae7edd1
--- /dev/null
+++ b/matlab/perfect-foresight-models/private/initialize_stacked_problem.m
@@ -0,0 +1,34 @@
+function [options, y0, yT, z, i_cols, i_cols_J1, i_cols_T, i_cols_j, i_cols_1, dynamicmodel] = initialize_stack_solve_algo_7(endogenousvariables, options, M)
+
+% Copyright (C) 2015 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 <http://www.gnu.org/licenses/>.
+
+if ~isfield(options.lmmcp,'lb')
+ [lb, ub, pfm.eq_index] = get_complementarity_conditions(M, options.ramsey_policy);
+ options.lmmcp.lb = repmat(lb, options.periods, 1);
+ options.lmmcp.ub = repmat(ub, options.periods, 1);
+end
+
+y0 = endogenousvariables(:,1);
+yT = endogenousvariables(:,options.periods+2);
+z = endogenousvariables(:,2:options.periods+1);
+illi = M.lead_lag_incidence';
+[i_cols, junk,i_cols_j] = find(illi(:));
+illi = illi(:,2:3);
+[i_cols_J1, junk,i_cols_1] = find(illi(:));
+i_cols_T = nonzeros(M.lead_lag_incidence(1:2,:)');
+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 85ecf6ac2a229a3c4b00b7a1d77afe239a081ff6..b52ecad172b4fc167e7dd48a00f7536532ece602 100644
--- a/matlab/perfect-foresight-models/private/simulation_core.m
+++ b/matlab/perfect-foresight-models/private/simulation_core.m
@@ -81,72 +81,12 @@ else
[oo_.endo_simul, oo_.deterministic_simulation] = ...
sim1_lbj(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
elseif options_.stack_solve_algo == 7
- periods = options_.periods;
- if ~isfield(options_.lmmcp,'lb')
- [lb,ub,pfm.eq_index] = get_complementarity_conditions(M_,options_.ramsey_policy);
- options_.lmmcp.lb = repmat(lb,periods,1);
- options_.lmmcp.ub = repmat(ub,periods,1);
- end
- y = oo_.endo_simul;
- y0 = y(:,1);
- yT = y(:,periods+2);
- z = y(:,2:periods+1);
- illi = M_.lead_lag_incidence';
- [i_cols,junk,i_cols_j] = find(illi(:));
- illi = illi(:,2:3);
- [i_cols_J1,junk,i_cols_1] = find(illi(:));
- i_cols_T = nonzeros(M_.lead_lag_incidence(1:2,:)');
- if options_.linear_approximation
- y_steady_state = oo_.steady_state;
- x_steady_state = transpose(oo_.exo_steady_state);
- ip = find(M_.lead_lag_incidence(1,:)');
- ic = find(M_.lead_lag_incidence(2,:)');
- in = find(M_.lead_lag_incidence(3,:)');
- % Evaluate the Jacobian of the dynamic model at the deterministic steady state.
- model_dynamic = str2func([M_.fname,'_dynamic']);
- [d1,jacobian] = model_dynamic(y_steady_state([ip; ic; in]), x_steady_state, M_.params, y_steady_state, 1);
- % Check that the dynamic model was evaluated at the steady state.
- if max(abs(d1))>1e-12
- error('Jacobian is not evaluated at the steady state!')
- end
- nyp = nnz(M_.lead_lag_incidence(1,:)) ;
- ny0 = nnz(M_.lead_lag_incidence(2,:)) ;
- nyf = nnz(M_.lead_lag_incidence(3,:)) ;
- nd = nyp+ny0+nyf; % size of y (first argument passed to the dynamic file).
- jexog = transpose(nd+(1:M_.exo_nbr));
- jendo = transpose(1:nd);
- z = bsxfun(@minus,z,y_steady_state);
- x = bsxfun(@minus,oo_.exo_simul,x_steady_state);
- [y,info] = dynare_solve(@linear_perfect_foresight_problem,z(:), options_, ...
- jacobian, y0-y_steady_state, yT-y_steady_state, ...
- x, M_.params, y_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_.NNZDerivatives(1),jendo,jexog);
- else
- [y,info] = dynare_solve(@perfect_foresight_problem,z(:),options_, ...
- 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_.NNZDerivatives(1));
- end
- if all(imag(y)<.1*options_.dynatol.f)
- if ~isreal(y)
- y = real(y);
- end
- else
- info = 1;
- end
if options_.linear_approximation
- oo_.endo_simul = [y0 bsxfun(@plus,reshape(y,M_.endo_nbr,periods),y_steady_state) yT];
- else
- oo_.endo_simul = [y0 reshape(y,M_.endo_nbr,periods) yT];
- end
- if info == 1
- oo_.deterministic_simulation.status = false;
+ [oo_.endo_simul, oo_.deterministic_simulation] = ...
+ solve_stacked_linear_problem(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
else
- oo_.deterministic_simulation.status = true;
+ [oo_.endo_simul, oo_.deterministic_simulation] = ...
+ solve_stacked_problem(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
end
end
end
diff --git a/matlab/perfect-foresight-models/solve_stacked_linear_problem.m b/matlab/perfect-foresight-models/solve_stacked_linear_problem.m
new file mode 100644
index 0000000000000000000000000000000000000000..b0f91b349780de44ac5f391fbc3cdef922f65a17
--- /dev/null
+++ b/matlab/perfect-foresight-models/solve_stacked_linear_problem.m
@@ -0,0 +1,65 @@
+function [endogenousvariables, info] = solve_stacked_linear_problem(endogenousvariables, exogenousvariables, steadystate_y, steadystate_x, M, options)
+
+% Copyright (C) 2015 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 <http://www.gnu.org/licenses/>.
+
+[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);
+
+ip = find(M_.lead_lag_incidence(1,:)');
+ic = find(M_.lead_lag_incidence(2,:)');
+in = find(M_.lead_lag_incidence(3,:)');
+
+% Evaluate the Jacobian of the dynamic model at the deterministic steady state.
+[d1,jacobian] = dynamicmodel(y_steady_state([ip; ic; in]), transpose(steadystate_x), M.params, steadystate_y, 1);
+
+% Check that the dynamic model was evaluated at the steady state.
+if max(abs(d1))>1e-12
+ error('Jacobian is not evaluated at the steady state!')
+end
+
+nyp = nnz(M_.lead_lag_incidence(1,:));
+ny0 = nnz(M_.lead_lag_incidence(2,:));
+nyf = nnz(M_.lead_lag_incidence(3,:));
+nd = nyp+ny0+nyf; % size of y (first argument passed to the dynamic file).
+jexog = transpose(nd+(1:M_.exo_nbr));
+jendo = transpose(1:nd);
+z = bsxfun(@minus, z, steadystate_y);
+x = bsxfun(@minus, exogenousvariables, steadystate_x);
+
+[y, check] = dynare_solve(@linear_perfect_foresight_problem,z(:), options, ...
+ 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, ...
+ M.NNZDerivatives(1), jendo, jexog);
+
+if all(imag(y)<.1*options.dynatol.x)
+ if ~isreal(y)
+ y = real(y);
+ end
+else
+ check = 1;
+end
+
+endogenousvariables = [y0 bsxfun(@plus,reshape(y,M.endo_nbr,options.periods), steadystate_y) yT];
+
+if check
+ info.status = false;
+else
+ info.status = true;
+end
\ No newline at end of file
diff --git a/matlab/perfect-foresight-models/solve_stacked_problem.m b/matlab/perfect-foresight-models/solve_stacked_problem.m
new file mode 100644
index 0000000000000000000000000000000000000000..7b340138c95bc42ad669c1e6b10629bbd580ed53
--- /dev/null
+++ b/matlab/perfect-foresight-models/solve_stacked_problem.m
@@ -0,0 +1,44 @@
+function [endogenousvariables, info] = solve_stacked_problem(endogenousvariables, exogenousvariables, steadystate, M, options);
+
+% Copyright (C) 2015 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 <http://www.gnu.org/licenses/>.
+
+[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);
+
+[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, ...
+ M.NNZDerivatives(1));
+
+if all(imag(y)<.1*options.dynatol.x)
+ if ~isreal(y)
+ y = real(y);
+ end
+else
+ check = 1;
+end
+
+endogenousvariables = [y0 reshape(y, M.endo_nbr, options.periods) yT];
+
+if check
+ info.status = false;
+else
+ info.status = true;
+end
\ No newline at end of file