extended_path.m 14.5 KB
Newer Older
1
2
function time_series = extended_path(initial_conditions,sample_size)
% Stochastic simulation of a non linear DSGE model using the Extended Path method (Fair and Taylor 1983). A time
3
4
% series of size T  is obtained by solving T perfect foresight models.
%
5
6
7
8
% INPUTS
%  o initial_conditions     [double]    m*nlags array, where m is the number of endogenous variables in the model and
%                                       nlags is the maximum number of lags.
%  o sample_size            [integer]   scalar, size of the sample to be simulated.
9
%
10
11
% OUTPUTS
%  o time_series            [double]    m*sample_size array, the simulations.
12
%
13
% ALGORITHM
14
%
15
16
% SPECIAL REQUIREMENTS

Sébastien Villemot's avatar
Sébastien Villemot committed
17
% Copyright (C) 2009-2013 Dynare Team
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
%
% 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/>.
global M_ options_ oo_
34

35
options_.verbosity = options_.ep.verbosity;
36
verbosity = options_.ep.verbosity+options_.ep.debug;
37

38
% Prepare a structure needed by the matlab implementation of the perfect foresight model solver
39
pfm = setup_stochastic_perfect_foresight_model_solver(M_,options_,oo_,'Tensor-Gaussian-Quadrature');
40

41
42
43
44
45
46
exo_nbr = M_.exo_nbr;
periods = options_.periods;
ep = options_.ep;
steady_state = oo_.steady_state;
dynatol = options_.dynatol;

47
48
49
50
51
52
% Set default initial conditions.
if isempty(initial_conditions)
    initial_conditions = oo_.steady_state;
end

% Set maximum number of iterations for the deterministic solver.
53
options_.simul.maxit = options_.ep.maxit;
54
55

% Set the number of periods for the perfect foresight model
56
periods = options_.ep.periods;
57
58
pfm.periods = options_.ep.periods;
pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny);
59

60
61
62
% keep a copy of pfm.i_upd
i_upd = pfm.i_upd;

63
64
65
% Set the algorithm for the perfect foresight solver
options_.stack_solve_algo = options_.ep.stack_solve_algo;

66
67
68
% Set check_stability flag
do_not_check_stability_flag = ~options_.ep.check_stability;

69
70
71
72
% Compute the first order reduced form if needed.
%
% REMARK. It is assumed that the user did run the same mod file with stoch_simul(order=1) and save
% all the globals in a mat file called linear_reduced_form.mat;
73

74
dr = struct();
75
if options_.ep.init
76
77
    options_.order = 1;
    [dr,Info,M_,options_,oo_] = resol(1,M_,options_,oo_);
78
79
80
end

% Do not use a minimal number of perdiods for the perfect foresight solver (with bytecode and blocks)
81
options_.minimal_solving_period = 100;%options_.ep.periods;
82
83
84
85
86
87
88
89

% Initialize the exogenous variables.
make_ex_;

% Initialize the endogenous variables.
make_y_;

% Initialize the output array.
90
time_series = zeros(M_.endo_nbr,sample_size);
91
92

% Set the covariance matrix of the structural innovations.
93
94
variances = diag(M_.Sigma_e);
positive_var_indx = find(variances>0);
95
96
effective_number_of_shocks = length(positive_var_indx);
stdd = sqrt(variances(positive_var_indx));
97
98
covariance_matrix = M_.Sigma_e(positive_var_indx,positive_var_indx);
covariance_matrix_upper_cholesky = chol(covariance_matrix);
99

100
% (re)Set exo_nbr
101
%exo_nbr = effective_number_of_shocks;
102

103
104
105
106
107
% Set seed.
if options_.ep.set_dynare_seed_to_default
    set_dynare_seed('default');
end

108
109
110
% Set bytecode flag
bytecode_flag = options_.ep.use_bytecode;

111
112
113
% Simulate shocks.
switch options_.ep.innovation_distribution
  case 'gaussian'
114
      oo_.ep.shocks = randn(sample_size,effective_number_of_shocks)*covariance_matrix_upper_cholesky;
115
116
117
  otherwise
    error(['extended_path:: ' options_.ep.innovation_distribution ' distribution for the structural innovations is not (yet) implemented!'])
end
118

119
120
121
% Initializes some variables.
t  = 0;

Stéphane Adjemian's avatar
Stéphane Adjemian committed
122
% Set waitbar (graphic or text  mode)
123
124
hh = dyn_waitbar(0,'Please wait. Extended Path simulations...');
set(hh,'Name','EP simulations.');
125

126
127
128
129
130
131
132
133
134
135
136
% hybrid correction
pfm.hybrid_order = options_.ep.stochastic.hybrid_order;
if pfm.hybrid_order
    oo_.dr = set_state_space(oo_.dr,M_,options_);
    options = options_;
    options.order = pfm.hybrid_order;
    pfm.dr = resol(0,M_,options,oo_);
else
    pfm.dr = [];
end

137
138
% Main loop.
while (t<sample_size)
139
    if ~mod(t,10)
140
        dyn_waitbar(t/sample_size,hh,'Please wait. Extended Path simulations...');
141
    end
142
143
144
145
146
    % Set period index.
    t = t+1;
    shocks = oo_.ep.shocks(t,:);
    % Put it in oo_.exo_simul (second line).
    oo_.exo_simul(2,positive_var_indx) = shocks;
147
148
149
150
151
152
153
154
155
156
157
158
    periods1 = periods;
    exo_simul_1 = zeros(periods1+2,exo_nbr);
    exo_simul_1(2,:) = oo_.exo_simul(2,:);
    pfm1 = pfm;
    info_convergence = 0;
    if ep.init% Compute first order solution (Perturbation)...
        ex = zeros(size(endo_simul_1,2),size(exo_simul_1,2));
        ex(1:size(exo_simul_1,1),:) = exo_simul_1;
        exo_simul_1 = ex;
        initial_path = simult_(initial_conditions,dr,exo_simul_1(2:end,:),1);
        endo_simul_1(:,1:end-1) = initial_path(:,1:end-1)*ep.init+endo_simul_1(:,1:end-1)*(1-ep.init);
    else
159
160
161
        if t==1
            endo_simul_1 = repmat(steady_state,1,periods1+2);
        end
162
163
164
    end
    % Solve a perfect foresight model.
    increase_periods = 0;
165
    % Keep a copy of endo_simul_1
166
167
168
    endo_simul = endo_simul_1;
    while 1
        if ~increase_periods
169
170
            if bytecode_flag && ~options_.ep.stochastic.order
                [flag,tmp] = bytecode('dynamic',endo_simul_1,exo_simul_1);
171
172
            else
                flag = 1;
173
            end
174
            if flag
Stéphane Adjemian's avatar
Stéphane Adjemian committed
175
176
                if options_.ep.stochastic.order == 0
                    [flag,tmp,err] = solve_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1);
177
                else
178
179
180
181
182
183
184
185
                    switch(options_.ep.stochastic.algo)
                        case 0
                        [flag,tmp] = ...
                            solve_stochastic_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1,options_.ep.stochastic.quadrature.nodes,options_.ep.stochastic.order);
                        case 1
                          [flag,tmp] = ...
                              solve_stochastic_perfect_foresight_model_1(endo_simul_1,exo_simul_1,pfm1,options_.ep.stochastic.quadrature.nodes,options_.ep.stochastic.order);
                    end
186
                end
187
188
            end
            info_convergence = ~flag;
189
        end
190
191
192
193
194
195
196
197
        if verbosity
            if info_convergence
                if t<10
                    disp(['Time:    ' int2str(t)  '. Convergence of the perfect foresight model solver!'])
                elseif t<100
                    disp(['Time:   ' int2str(t)  '. Convergence of the perfect foresight model solver!'])
                elseif t<1000
                    disp(['Time:  ' int2str(t)  '. Convergence of the perfect foresight model solver!'])
198
                else
199
                    disp(['Time: ' int2str(t)  '. Convergence of the perfect foresight model solver!'])
200
                end
201
202
203
204
205
206
207
208
209
            else
                if t<10
                    disp(['Time:    ' int2str(t)  '. No convergence of the perfect foresight model solver!'])
                elseif t<100
                    disp(['Time:   ' int2str(t)  '. No convergence of the perfect foresight model solver!'])
                elseif t<1000
                    disp(['Time:  ' int2str(t)  '. No convergence of the perfect foresight model solver!'])
                else
                    disp(['Time: ' int2str(t)  '. No convergence of the perfect foresight model solver!'])
210
                end
211
            end
212
213
214
215
216
217
218
219
220
221
222
223
224
        end
        if do_not_check_stability_flag
            % Exit from the while loop.
            endo_simul_1 = tmp;
            break
        else
            % Test if periods is big enough.
            % Increase the number of periods.
            periods1 = periods1 + ep.step;
            pfm1.periods = periods1;
            pfm1.i_upd = pfm1.ny+(1:pfm1.periods*pfm1.ny);
            % Increment the counter.
            increase_periods = increase_periods + 1;
225
            if verbosity
226
227
228
229
230
231
                if t<10
                    disp(['Time:    ' int2str(t)  '. I increase the number of periods to ' int2str(periods1) '.'])
                elseif t<100
                    disp(['Time:   ' int2str(t) '. I increase the number of periods to ' int2str(periods1) '.'])
                elseif t<1000
                    disp(['Time:  ' int2str(t)  '. I increase the number of periods to ' int2str(periods1) '.'])
232
                else
233
                    disp(['Time: ' int2str(t)  '. I increase the number of periods to ' int2str(periods1) '.'])
234
235
                end
            end
236
237
238
239
            if info_convergence
                % If the previous call to the perfect foresight model solver exited
                % announcing that the routine converged, adapt the size of endo_simul_1
                % and exo_simul_1.
Stéphane Adjemian's avatar
Stéphane Adjemian committed
240
241
                endo_simul_1 = [ tmp , repmat(steady_state,1,ep.step) ];
                exo_simul_1  = [ exo_simul_1 ; zeros(ep.step,exo_nbr)];
242
                tmp_old = tmp;
243
            else
244
245
246
247
248
249
                % If the previous call to the perfect foresight model solver exited
                % announcing that the routine did not converge, then tmp=1... Maybe
                % should change that, because in some circonstances it may usefull
                % to know where the routine did stop, even if convergence was not
                % achieved.
                endo_simul_1 = [ endo_simul_1 , repmat(steady_state,1,ep.step) ];
Stéphane Adjemian's avatar
Stéphane Adjemian committed
250
                exo_simul_1  = [ exo_simul_1 ; zeros(ep.step,exo_nbr)];
251
252
            end
            % Solve the perfect foresight model with an increased number of periods.
253
254
            if bytecode_flag && ~options_.ep.stochastic.order
                [flag,tmp] = bytecode('dynamic',endo_simul_1,exo_simul_1);
255
256
257
258
            else
                flag = 1;
            end
            if flag
Stéphane Adjemian's avatar
Stéphane Adjemian committed
259
260
                if options_.ep.stochastic.order == 0
                    [flag,tmp,err] = solve_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1);
261
                else
Stéphane Adjemian's avatar
Stéphane Adjemian committed
262
                    [flag,tmp] = solve_stochastic_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1,options_.ep.stochastic.nodes,options_.ep.stochastic.order);
263
                end
264
265
266
267
268
269
270
            end
            info_convergence = ~flag;
            if info_convergence
                % If the solver achieved convergence, check that simulated paths did not
                % change during the first periods.
                % Compute the maximum deviation between old path and new path over the
                % first periods
Stéphane Adjemian's avatar
Stéphane Adjemian committed
271
                delta = max(max(abs(tmp(:,2)-tmp_old(:,2))));
272
273
274
275
276
277
278
279
280
281
282
                if delta < dynatol.x
                    % If the maximum deviation is close enough to zero, reset the number
                    % of periods to ep.periods
                    periods1 = ep.periods;
                    pfm1.periods = periods1;
                    pfm1.i_upd = pfm1.ny+(1:pfm1.periods*pfm1.ny);
                    % Cut exo_simul_1 and endo_simul_1 consistently with the resetted
                    % number of periods and exit from the while loop.
                    exo_simul_1 = exo_simul_1(1:(periods1+2),:);
                    endo_simul_1 = endo_simul_1(:,1:(periods1+2));
                    break
283
                end
284
285
286
287
288
289
290
291
292
293
294
295
296
297
            else
                % The solver did not converge... Try to solve the model again with a bigger
                % number of periods, except if the number of periods has been increased more
                % than 10 times.
                if increase_periods==10;
                    if verbosity
                        if t<10
                            disp(['Time:    ' int2str(t)  '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
                        elseif t<100
                            disp(['Time:   ' int2str(t)  '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
                        elseif t<1000
                            disp(['Time:  ' int2str(t)  '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
                        else
                            disp(['Time: ' int2str(t)  '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
298
                        end
299
                    end
300
301
                    % Exit from the while loop.
                    break
302
                end
303
304
305
306
            end% if info_convergence
        end
    end% while
    if ~info_convergence% If exited from the while loop without achieving convergence, use an homotopic approach
307
308
309
310
311
312
313
        if ~do_not_check_stability_flag
            periods1 = ep.periods;
            pfm1.periods = periods1;
            pfm1.i_upd = i_upd;
            exo_simul_1 = exo_simul_1(1:(periods1+2),:);
            endo_simul_1 = endo_simul_1(:,1:(periods1+2));
        end
Stéphane Adjemian's avatar
Stéphane Adjemian committed
314
        [INFO,tmp] = homotopic_steps(endo_simul,exo_simul_1,.5,.01,pfm1);
315
316
317
318
319
320
        if isstruct(INFO)
            info_convergence = INFO.convergence;
        else
            info_convergence = 0;
        end
        if ~info_convergence
Stéphane Adjemian's avatar
Stéphane Adjemian committed
321
            [INFO,tmp] = homotopic_steps(endo_simul,exo_simul_1,0,.01,pfm1);
322
323
324
325
326
            if isstruct(INFO)
                info_convergence = INFO.convergence;
            else
                info_convergence = 0;
            end
327
328
329
            if ~info_convergence
                disp('Homotopy:: No convergence of the perfect foresight model solver!')
                error('I am not able to simulate this model!');
330
            else
331
332
                endo_simul_1 = tmp;
                if verbosity && info_convergence
333
334
335
                    disp('Homotopy:: Convergence of the perfect foresight model solver!')
                end
            end
336
337
338
339
340
341
        else
            info_convergence = 1;
            endo_simul_1 = tmp;
            if verbosity && info_convergence
                disp('Homotopy:: Convergence of the perfect foresight model solver!')
            end
342
343
        end
    end
344
    % Save results of the perfect foresight model solver.
345
346
347
348
    time_series(:,t) = endo_simul_1(:,2);
    endo_simul_1(:,1:end-1) = endo_simul_1(:,2:end);
    endo_simul_1(:,1) = time_series(:,t);
    endo_simul_1(:,end) = oo_.steady_state;
349
end% (while) loop over t
350

351
dyn_waitbar_close(hh);
352

353
oo_.endo_simul = oo_.steady_state;