perfect_foresight_solver_core.m 6.04 KB
Newer Older
1
function [oo_, maxerror] = perfect_foresight_solver_core(M_, options_, oo_)
2

3 4
% Core function calling solvers for perfect foresight model
%
5
% INPUTS
6 7 8 9
% - M_                  [struct] contains a description of the model.
% - options_            [struct] contains various options.
% - oo_                 [struct] contains results
%
10
% OUTPUTS
11 12
% - oo_                 [struct] contains results
% - maxerror            [double] contains the maximum absolute error
13

14
% Copyright (C) 2015-2019 Dynare Team
15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
%
% 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/>.

Stéphane Adjemian's avatar
Stéphane Adjemian committed
31 32 33 34 35
if options_.lmmcp.status
    options_.stack_solve_algo=7;
    options_.solve_algo = 10;
end

36 37
periods = options_.periods;

38
if 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.')
40 41
end

42 43
if options_.linear && (isequal(options_.stack_solve_algo, 0) || isequal(options_.stack_solve_algo, 7))
    options_.linear_approximation = true;
44
end
45

46 47 48
if options_.block
    if options_.bytecode
        try
49
            [info, tmp] = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state,1, periods+2), periods);
50
        catch
51
            info = 1;
52 53
        end
        if info
54
            oo_.deterministic_simulation.status = false;
55
        else
56 57 58 59 60
            oo_.endo_simul = tmp;
            oo_.deterministic_simulation.status = true;
        end
        if options_.no_homotopy
            mexErrCheck('bytecode', info);
61 62
        end
    else
63
        oo_ = feval([M_.fname '.dynamic'], options_, M_, oo_);
64 65 66 67
    end
else
    if options_.bytecode
        try
68
            [info, tmp] = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state, 1, periods+2), periods);
69
        catch
70
            info = 1;
71 72
        end
        if info
73
            oo_.deterministic_simulation.status = false;
74
        else
75 76 77 78 79 80
            oo_.endo_simul = tmp;
            oo_.deterministic_simulation.status = true;
        end
        if options_.no_homotopy
            mexErrCheck('bytecode', info);
        end
81 82
    else
        if M_.maximum_endo_lead == 0 % Purely backward model
83 84
            [oo_.endo_simul, oo_.deterministic_simulation] = ...
                sim1_purely_backward(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
85
        elseif M_.maximum_endo_lag == 0 % Purely forward model
86 87
        [oo_.endo_simul, oo_.deterministic_simulation] = ...
            sim1_purely_forward(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
88
        else % General case
Stéphane Adjemian's avatar
Stéphane Adjemian committed
89 90
            switch options_.stack_solve_algo
              case 0
91
                if options_.linear_approximation
92 93
                    [oo_.endo_simul, oo_.deterministic_simulation] = ...
                        sim1_linear(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, oo_.exo_steady_state, M_, options_);
94
                else
95 96
                    [oo_.endo_simul, oo_.deterministic_simulation] = ...
                        sim1(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
97
                end
Stéphane Adjemian's avatar
Stéphane Adjemian committed
98
              case 6
Stéphane Adjemian's avatar
Stéphane Adjemian committed
99
                if options_.linear_approximation
Stéphane Adjemian's avatar
Stéphane Adjemian committed
100 101
                    error('Invalid value of stack_solve_algo option!')
                end
102 103
                [oo_.endo_simul, oo_.deterministic_simulation] = ...
                    sim1_lbj(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
Stéphane Adjemian's avatar
Stéphane Adjemian committed
104
              case 7
105
                if options_.linear_approximation
Stéphane Adjemian's avatar
Stéphane Adjemian committed
106 107 108
                    if isequal(options_.solve_algo, 10)
                        warning('It would be more efficient to set option solve_algo equal to 0!')
                    end
Stéphane Adjemian's avatar
Stéphane Adjemian committed
109 110
                    [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_);
111
                else
Stéphane Adjemian's avatar
Stéphane Adjemian committed
112 113
                    [oo_.endo_simul, oo_.deterministic_simulation] = ...
                        solve_stacked_problem(oo_.endo_simul, oo_.exo_simul, oo_.steady_state, M_, options_);
114
                end
Stéphane Adjemian's avatar
Stéphane Adjemian committed
115 116
              otherwise
                error('Invalid value of stack_solve_algo option!')
117 118 119 120 121
            end
        end
    end
end

122
if nargout>1
123
    y0 = oo_.endo_simul(:,1);
124 125 126 127 128 129 130 131 132 133 134 135 136
    yT = oo_.endo_simul(:,periods+2);
    yy  = oo_.endo_simul(:,2:periods+1);
    illi = M_.lead_lag_incidence';
    [i_cols,~,i_cols_j] = find(illi(:));
    illi = illi(:,2:3);
    [i_cols_J1,~,i_cols_1] = find(illi(:));
    i_cols_T = nonzeros(M_.lead_lag_incidence(1:2,:)');
    if periods==1
        i_cols_0 = nonzeros(M_.lead_lag_incidence(2,:)');
        i_cols_J0 = find(M_.lead_lag_incidence(2,:)');
    else
        i_cols_0 = [];
        i_cols_J0 = [];
137
    end
138 139 140 141
    if options_.block && ~options_.bytecode
        maxerror = oo_.deterministic_simulation.error;
    else
        if options_.bytecode
142
            [chck, residuals, ~]= bytecode('dynamic','evaluate', oo_.endo_simul, oo_.exo_simul, M_.params, oo_.steady_state, 1);
143
        else
144
            residuals = perfect_foresight_problem(yy(:),str2func([M_.fname '.dynamic']), y0, yT, ...
145
                                                  oo_.exo_simul,M_.params,oo_.steady_state, ...
146
                                                  M_.maximum_lag, periods,M_.endo_nbr,i_cols, ...
147
                                                  i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0);
148 149 150
        end
        maxerror = max(max(abs(residuals)));
    end
151
end