1function [flag,endo_simul,err] = solve_perfect_foresight_model(endo_simul,exo_simul,pfm)
2
3% Copyright (C) 2012-2017 Dynare Team
4%
5% This file is part of Dynare.
6%
7% Dynare is free software: you can redistribute it and/or modify
8% it under the terms of the GNU General Public License as published by
9% the Free Software Foundation, either version 3 of the License, or
10% (at your option) any later version.
11%
12% Dynare is distributed in the hope that it will be useful,
13% but WITHOUT ANY WARRANTY; without even the implied warranty of
14% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15% GNU General Public License for more details.
16%
17% You should have received a copy of the GNU General Public License
18% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
19
20flag = 0;
21err = 0;
22stop = 0;
23nan_flag = 0;
24
25model_dynamic = pfm.dynamic_model;
26
27Y = endo_simul(:);
28
29if pfm.verbose
30    disp (['-----------------------------------------------------']) ;
31    disp (['MODEL SIMULATION :']) ;
32    fprintf('\n') ;
33end
34
35if pfm.use_bytecode
36    [flag, endo_simul]=bytecode(Y, exo_simul, pfm.params);
37    return
38end
39
40z = Y(find(pfm.lead_lag_incidence'));
41[d1,jacobian] = model_dynamic(z,exo_simul,pfm.params,pfm.steady_state,2);
42
43% Initialization of the jacobian of the stacked model.
44A = sparse([],[],[],pfm.periods*pfm.ny,pfm.periods*pfm.ny,pfm.periods*nnz(jacobian));
45
46% Initialization of the Newton residuals.
47res = zeros(pfm.periods*pfm.ny,1);
48
49h1 = clock;
50
51% Newton loop.
52for iter = 1:pfm.maxit_
53    h2 = clock;
54    i_rows = 1:pfm.ny;
55    i_cols = find(pfm.lead_lag_incidence');
56    i_cols_A = i_cols;
57    % Fill the jacobian of the stacked model.
58    for it = 2:(pfm.periods+1)
59        [d1,jacobian] = model_dynamic(Y(i_cols),exo_simul,pfm.params,pfm.steady_state,it);
60        if it == 2
61            A(i_rows,pfm.i_cols_A1) = jacobian(:,pfm.i_cols_1);
62        elseif it == pfm.periods+1
63            A(i_rows,i_cols_A(pfm.i_cols_T)) = jacobian(:,pfm.i_cols_T);
64        else
65            A(i_rows,i_cols_A) = jacobian(:,pfm.i_cols_j);
66        end
67        res(i_rows) = d1;
68        i_rows = i_rows + pfm.ny;
69        i_cols = i_cols + pfm.ny;
70        if it > 2
71            i_cols_A = i_cols_A + pfm.ny;
72        end
73    end
74    % Stop if Newton residuals are zero.
75    err = max(abs(res));
76    if err < pfm.tolerance
77        stop = 1 ;
78        if pfm.verbose
79            fprintf('\n') ;
80            disp([' Total time of simulation        :' num2str(etime(clock,h1))]) ;
81            fprintf('\n') ;
82            disp([' Convergency obtained.']) ;
83            fprintf('\n') ;
84        end
85        flag = 0;% Convergency obtained.
86        endo_simul = reshape(Y,pfm.ny,pfm.periods+2);
87        break
88    end
89    % Compute the Newton step.
90    dy = -A\res;
91    if any(isnan(dy))
92        nan_flag = 1;
93        break
94    end
95    % Update the endogenous variables paths.
96    Y(pfm.i_upd) =   Y(pfm.i_upd) + dy;
97end
98
99if ~stop
100    if pfm.verbose
101        fprintf('\n') ;
102        disp(['     Total time of simulation        :' num2str(etime(clock,h1))]) ;
103        fprintf('\n') ;
104        disp(['WARNING : maximum number of iterations is reached (modify options_.simul.maxit).']) ;
105        fprintf('\n') ;
106    end
107    flag = 1;% more iterations are needed.
108    endo_simul = 1;
109end
110if nan_flag
111    if pfm.verbose
112        fprintf('\n') ;
113        disp(['     Total time of simulation        :' num2str(etime(clock,h1))]) ;
114        fprintf('\n') ;
115        disp(['WARNING : NaNs!']) ;
116        fprintf('\n') ;
117    end
118    flag = 1;
119    endo_simul = 1;
120end
121if pfm.verbose
122    disp (['-----------------------------------------------------']) ;
123end