1function [oo_, maxerror] = perfect_foresight_solver_core(M_, options_, oo_)
2
3% Core function calling solvers for perfect foresight model
4%
5% INPUTS
6% - M_                  [struct] contains a description of the model.
7% - options_            [struct] contains various options.
8% - oo_                 [struct] contains results
9%
10% OUTPUTS
11% - oo_                 [struct] contains results
12% - maxerror            [double] contains the maximum absolute error
13
14% Copyright (C) 2015-2019 Dynare Team
15%
16% This file is part of Dynare.
17%
18% Dynare is free software: you can redistribute it and/or modify
19% it under the terms of the GNU General Public License as published by
20% the Free Software Foundation, either version 3 of the License, or
21% (at your option) any later version.
22%
23% Dynare is distributed in the hope that it will be useful,
24% but WITHOUT ANY WARRANTY; without even the implied warranty of
25% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
26% GNU General Public License for more details.
27%
28% You should have received a copy of the GNU General Public License
29% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
30
31if options_.lmmcp.status
32    options_.stack_solve_algo=7;
33    options_.solve_algo = 10;
34end
35
36periods = options_.periods;
37
38if options_.linear_approximation && ~(isequal(options_.stack_solve_algo,0) || isequal(options_.stack_solve_algo,7))
39    error('perfect_foresight_solver: Option linear_approximation is only available with option stack_solve_algo equal to 0 or 7.')
40end
41
42if options_.endogenous_terminal_period && options_.stack_solve_algo ~= 0
43    error('perfect_foresight_solver: option endogenous_terminal_period is only available with option stack_solve_algo equal to 0')
44end
45
46if options_.linear && (isequal(options_.stack_solve_algo, 0) || isequal(options_.stack_solve_algo, 7))
47    options_.linear_approximation = true;
48end
49
50if options_.block
51    if options_.bytecode
52        try
53            [info, tmp] = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state,1, periods+2), periods);
54        catch
55            info = 1;
56        end
57        if info
58            oo_.deterministic_simulation.status = false;
59        else
60            oo_.endo_simul = tmp;
61            oo_.deterministic_simulation.status = true;
62        end
63        if options_.no_homotopy
64            mexErrCheck('bytecode', info);
65        end
66    else
67        oo_ = feval([M_.fname '.dynamic'], options_, M_, oo_);
68    end
69else
70    if options_.bytecode
71        try
72            [info, tmp] = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state, 1, periods+2), periods);
73        catch
74            info = 1;
75        end
76        if info
77            oo_.deterministic_simulation.status = false;
78        else
79            oo_.endo_simul = tmp;
80            oo_.deterministic_simulation.status = true;
81        end
82        if options_.no_homotopy
83            mexErrCheck('bytecode', info);
84        end
85    else
86        if M_.maximum_endo_lead == 0 && ~options_.lmmcp.status % Purely backward model
87            [oo_.endo_simul, oo_.deterministic_simulation] = ...
88                sim1_purely_backward(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
89        elseif M_.maximum_endo_lag == 0 && ~options_.lmmcp.status % Purely forward model
90        [oo_.endo_simul, oo_.deterministic_simulation] = ...
91            sim1_purely_forward(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
92        else % General case
93            switch options_.stack_solve_algo
94              case 0
95                if options_.linear_approximation
96                    [oo_.endo_simul, oo_.deterministic_simulation] = ...
97                        sim1_linear(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
98                else
99                    [oo_.endo_simul, oo_.deterministic_simulation] = ...
100                        sim1(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
101                end
102              case 6
103                if options_.linear_approximation
104                    error('Invalid value of stack_solve_algo option!')
105                end
106                [oo_.endo_simul, oo_.deterministic_simulation] = ...
107                    sim1_lbj(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
108              case 7
109                if options_.linear_approximation
110                    if isequal(options_.solve_algo, 10)
111                        warning('It would be more efficient to set option solve_algo equal to 0!')
112                    end
113                    [oo_.endo_simul, oo_.deterministic_simulation] = ...
114                        solve_stacked_linear_problem(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
115                else
116                    [oo_.endo_simul, oo_.deterministic_simulation] = ...
117                        solve_stacked_problem(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
118                end
119              otherwise
120                error('Invalid value of stack_solve_algo option!')
121            end
122        end
123    end
124end
125
126if nargout>1
127    if options_.block && ~options_.bytecode
128        maxerror = oo_.deterministic_simulation.error;
129    else
130        if options_.bytecode
131            [~, residuals]= bytecode('dynamic','evaluate', oo_.endo_simul, oo_.exo_simul, M_.params, oo_.steady_state, 1);
132        else
133            if M_.maximum_lag > 0
134                y0 = oo_.endo_simul(:, M_.maximum_lag);
135            else
136                y0 = NaN(ny, 1);
137            end
138            if M_.maximum_lead > 0
139                yT = oo_.endo_simul(:, M_.maximum_lag+periods+1);
140            else
141                yT = NaN(ny, 1);
142            end
143            yy = oo_.endo_simul(:,M_.maximum_lag+(1:periods));
144
145            residuals = perfect_foresight_problem(yy(:), y0, yT, oo_.exo_simul, M_.params, oo_.steady_state, periods, M_, options_);
146        end
147        maxerror = max(max(abs(residuals)));
148    end
149end
150