1function perfect_foresight_solver()
2% Computes deterministic simulations
3%
4% INPUTS
5%   None
6%
7% OUTPUTS
8%   none
9%
10% ALGORITHM
11%
12% SPECIAL REQUIREMENTS
13%   none
14
15% Copyright (C) 1996-2020 Dynare Team
16%
17% This file is part of Dynare.
18%
19% Dynare is free software: you can redistribute it and/or modify
20% it under the terms of the GNU General Public License as published by
21% the Free Software Foundation, either version 3 of the License, or
22% (at your option) any later version.
23%
24% Dynare is distributed in the hope that it will be useful,
25% but WITHOUT ANY WARRANTY; without even the implied warranty of
26% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
27% GNU General Public License for more details.
28%
29% You should have received a copy of the GNU General Public License
30% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
31
32global M_ options_ oo_
33
34check_input_arguments(options_, M_, oo_);
35
36if isempty(options_.scalv) || options_.scalv == 0
37    options_.scalv = oo_.steady_state;
38end
39
40periods = options_.periods;
41
42options_.scalv= 1;
43
44if options_.debug
45    model_static = str2func([M_.fname,'.static']);
46    for ii=1:size(oo_.exo_simul,1)
47        [residual(:,ii)] = model_static(oo_.steady_state, oo_.exo_simul(ii,:),M_.params);
48    end
49    problematic_periods=find(any(isinf(residual)) | any(isnan(residual)))-M_.maximum_endo_lag;
50    if ~isempty(problematic_periods)
51        period_string=num2str(problematic_periods(1));
52        for ii=2:length(problematic_periods)
53            period_string=[period_string, ', ', num2str(problematic_periods(ii))];
54        end
55        fprintf('\n\nWARNING: Value for the exogenous variable(s) in period(s) %s inconsistent with the static model.\n',period_string);
56        fprintf('WARNING: Check for division by 0.\n')
57    end
58end
59
60initperiods = 1:M_.maximum_lag;
61lastperiods = (M_.maximum_lag+periods+1):(M_.maximum_lag+periods+M_.maximum_lead);
62
63oo_ = perfect_foresight_solver_core(M_,options_,oo_);
64
65% If simulation failed try homotopy.
66if ~oo_.deterministic_simulation.status && ~options_.no_homotopy
67
68    if ~options_.noprint
69        fprintf('\nSimulation of the perfect foresight model failed!')
70        fprintf('Switching to a homotopy method...\n')
71    end
72
73    if ~M_.maximum_lag
74        disp('Homotopy not implemented for purely forward models!')
75        disp('Failed to solve the model!')
76        disp('Return with empty oo_.endo_simul.')
77        oo_.endo_simul = [];
78        return
79    end
80    if ~M_.maximum_lead
81        disp('Homotopy not implemented for purely backward models!')
82        disp('Failed to solve the model!')
83        disp('Return with empty oo_.endo_simul.')
84        oo_.endo_simul = [];
85        return
86    end
87
88    % Disable warnings if homotopy
89    warning_old_state = warning;
90    warning off all
91    % Do not print anything
92    oldverbositylevel = options_.verbosity;
93    options_.verbosity = 0;
94
95    % Set initial paths for the endogenous and exogenous variables.
96    endoinit = repmat(oo_.steady_state, 1,M_.maximum_lag+periods+M_.maximum_lead);
97    exoinit = repmat(oo_.exo_steady_state',M_.maximum_lag+periods+M_.maximum_lead,1);
98
99    % Copy the current paths for the exogenous and endogenous variables.
100    exosim = oo_.exo_simul;
101    endosim = oo_.endo_simul;
102
103    current_weight = 0;    % Current weight of target point in convex combination.
104    step = .5;             % Set default step size.
105    success_counter = 0;
106    iteration = 0;
107
108    if ~options_.noprint
109        fprintf('Iter. \t | Lambda \t | status \t | Max. residual\n')
110        fprintf('++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n')
111    end
112    while (step > options_.dynatol.x)
113
114        if ~isequal(step,1)
115            options_.verbosity = 0;
116        end
117
118        iteration = iteration+1;
119        new_weight = current_weight + step; % Try this weight, and see if it succeeds
120
121        if new_weight >= 1
122            new_weight = 1; % Don't go beyond target point
123            step = new_weight - current_weight;
124        end
125
126        % Compute convex combination for exo path and initial/terminal endo conditions
127        % But take care of not overwriting the computed part of oo_.endo_simul
128        oo_.exo_simul = exosim*new_weight + exoinit*(1-new_weight);
129        oo_.endo_simul(:,[initperiods, lastperiods]) = new_weight*endosim(:,[initperiods, lastperiods])+(1-new_weight)*endoinit(:,[initperiods, lastperiods]);
130
131        % Detect Nans or complex numbers in the solution.
132        path_with_nans = any(any(isnan(oo_.endo_simul)));
133        path_with_cplx = any(any(~isreal(oo_.endo_simul)));
134
135        if isequal(iteration, 1)
136            % First iteration, same initial guess as in the first call to perfect_foresight_solver_core routine.
137            oo_.endo_simul(:,M_.maximum_lag+1:end-M_.maximum_lead) = endoinit(:,1:periods);
138        elseif path_with_nans || path_with_cplx
139            % If solver failed with NaNs or complex number, use previous solution as an initial guess.
140            oo_.endo_simul(:,M_.maximum_lag+1:end-M_.maximum_lead) = saved_endo_simul(:,1+M_.maximum_lag:end-M_.maximum_lead);
141        end
142
143        % Make a copy of the paths.
144        saved_endo_simul = oo_.endo_simul;
145
146        % Solve for the paths of the endogenous variables.
147        [oo_,me] = perfect_foresight_solver_core(M_,options_,oo_);
148
149        if oo_.deterministic_simulation.status == 1
150            current_weight = new_weight;
151            if current_weight >= 1
152                if ~options_.noprint
153                    fprintf('%i \t | %1.5f \t | %s \t | %e\n', iteration, new_weight, 'succeeded', me)
154                end
155                break
156            end
157            success_counter = success_counter + 1;
158            if success_counter >= 3
159                success_counter = 0;
160                step = step * 2;
161            end
162            if ~options_.noprint
163                fprintf('%i \t | %1.5f \t | %s \t | %e\n', iteration, new_weight, 'succeeded', me)
164            end
165        else
166            % If solver failed, then go back.
167            oo_.endo_simul = saved_endo_simul;
168            success_counter = 0;
169            step = step / 2;
170            if ~options_.noprint
171                if isreal(me)
172                    fprintf('%i \t | %1.5f \t | %s \t | %e\n', iteration, new_weight, 'failed', me)
173                else
174                    fprintf('%i \t | %1.5f \t | %s \t | %s\n', iteration, new_weight, 'failed', 'Complex')
175                end
176            end
177        end
178    end
179    if ~options_.noprint
180        fprintf('++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n\n')
181    end
182    options_.verbosity = oldverbositylevel;
183    warning(warning_old_state);
184end
185
186
187if ~isreal(oo_.endo_simul(:)) % cannot happen with bytecode or the perfect_foresight_problem DLL
188    ny = size(oo_.endo_simul, 1)
189    if M_.maximum_lag > 0
190        y0 = real(oo_.endo_simul(:, M_.maximum_lag));
191    else
192        y0 = NaN(ny, 1);
193    end
194    if M_.maximum_lead > 0
195        yT = real(oo_.endo_simul(:, M_.maximum_lag+periods+1));
196    else
197        yT = NaN(ny, 1);
198    end
199    yy = real(oo_.endo_simul(:,M_.maximum_lag+(1:periods)));
200
201    residuals = perfect_foresight_problem(yy(:), y0, yT, oo_.exo_simul, M_.params, oo_.steady_state, periods, M_, options_);
202
203    if max(abs(residuals))< options_.dynatol.f
204        oo_.deterministic_simulation.status = 1;
205        oo_.endo_simul=real(oo_.endo_simul);
206    else
207        oo_.deterministic_simulation.status = 0;
208        disp('Simulation terminated with imaginary parts in the residuals or endogenous variables.')
209    end
210end
211
212if oo_.deterministic_simulation.status == 1
213    if ~options_.noprint
214        fprintf('Perfect foresight solution found.\n\n')
215    end
216else
217    fprintf('Failed to solve perfect foresight model\n\n')
218end
219
220dyn2vec(M_, oo_, options_);
221
222if ~isdates(options_.initial_period) && isnan(options_.initial_period)
223    initial_period = dates(1,1);
224else
225    initial_period = options_.initial_period;
226end
227
228ts = dseries(transpose(oo_.endo_simul), initial_period, M_.endo_names);
229assignin('base', 'Simulated_time_series', ts);
230if oo_.deterministic_simulation.status
231    oo_.gui.ran_perfect_foresight = true;
232end
233