solve_stochastic_perfect_foresight_model.m 8.66 KB
Newer Older
1
function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model(endo_simul,exo_simul,pfm,nnodes,order)
2

Sébastien Villemot's avatar
Sébastien Villemot committed
3
% Copyright (C) 2012-2013 Dynare Team
Sébastien Villemot's avatar
Sébastien Villemot committed
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
%
% 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/>.

20
21
22
23
    flag = 0;
    err = 0;
    stop = 0;

24
25
26
27
28
29
30
31
32
33
34
35
    params = pfm.params;
    steady_state = pfm.steady_state;
    ny = pfm.ny;
    periods = pfm.periods;
    dynamic_model = pfm.dynamic_model;
    lead_lag_incidence = pfm.lead_lag_incidence;
    nyp = pfm.nyp;
    nyf = pfm.nyf;
    i_cols_1 = pfm.i_cols_1;
    i_cols_A1 = pfm.i_cols_A1;
    i_cols_j = pfm.i_cols_j;
    i_cols_T = nonzeros(lead_lag_incidence(1:2,:)');
36

37
38
39
    maxit = pfm.maxit_;
    tolerance = pfm.tolerance;
    verbose = pfm.verbose;
40

41
42
43
    number_of_shocks = size(exo_simul,2);

    [nodes,weights] = gauss_hermite_weights_and_nodes(nnodes);
44

45
    if number_of_shocks>1
46
        nodes = repmat(nodes,1,number_of_shocks)*chol(pfm.Sigma);
47
        % to be fixed for Sigma ~= I
48
        for i=1:number_of_shocks
49
            rr(i) = {nodes(:,i)};
50
51
52
53
            ww(i) = {weights};
        end
        nodes = cartesian_product_of_sets(rr{:});
        weights = prod(cartesian_product_of_sets(ww{:}),2);
54
        nnodes = nnodes^number_of_shocks;
55
    else
56
        nodes = nodes*sqrt(pfm.Sigma);
57
58
    end

59
    innovations = zeros(periods+2,number_of_shocks);
60

61
    if verbose
62
63
64
        disp ([' -----------------------------------------------------']);
        disp (['MODEL SIMULATION :']);
        fprintf('\n');
65
66
    end

67
68
    z = endo_simul(find(lead_lag_incidence'));
    [d1,jacobian] = dynamic_model(z,exo_simul,params,steady_state,2);
69

70
71
72
73
74
75
76
77
    % Each column of Y represents a different world
    % The upper right cells are unused
    % The first row block is ny x 1
    % The second row block is ny x nnodes
    % The third row block is ny x nnodes^2
    % and so on until size ny x nnodes^order
    world_nbr = nnodes^order;
    Y = repmat(endo_simul(:),1,world_nbr);
78

79
80
81
82
    % The columns of A map the elements of Y such that
    % each block of Y with ny rows are unfolded column wise
    dimension = ny*(sum(nnodes.^(0:order-1),2)+(periods-order)*world_nbr);
    if order == 0
83
84
        i_upd_r = (1:ny*periods);
        i_upd_y = i_upd_r + ny;
85
    else
86
87
88
89
        i_upd_r = zeros(dimension,1);
        i_upd_y = i_upd_r;
        i_upd_r(1:ny) = (1:ny);
        i_upd_y(1:ny) = ny+(1:ny);
90
        i1 = ny+1;
91
        i2 = 2*ny;
92
93
        n1 = ny+1;
        n2 = 2*ny;
94
95
96
        for i=2:periods
            k = n1:n2;
            for j=1:nnodes^min(i-1,order)
97
98
                i_upd_r(i1:i2) = (n1:n2)+(j-1)*ny*periods;
                i_upd_y(i1:i2) = (n1:n2)+ny+(j-1)*ny*(periods+2);
99
                i1 = i2+1;
100
                i2 = i2+ny;
101
            end
102
103
            n1 = n2+1;
            n2 = n2+ny;
104
105
        end
    end
106
107
    icA = [find(lead_lag_incidence(1,:)) find(lead_lag_incidence(2,:))+world_nbr*ny ...
           find(lead_lag_incidence(3,:))+2*world_nbr*ny]';
108
    h1 = clock;
109
    for iter = 1:maxit
110
        h2 = clock;
111
112
        A1 = sparse([],[],[],ny*(sum(nnodes.^(0:order-1),2)+1),dimension,(order+1)*world_nbr*nnz(jacobian));
        res = zeros(ny,periods,world_nbr);
113
114
115
        i_rows = 1:ny;
        i_cols = find(lead_lag_incidence');
        i_cols_p = i_cols(1:nyp);
116
117
        i_cols_s = i_cols(nyp+(1:ny));
        i_cols_f = i_cols(nyp+ny+(1:nyf));
118
        i_cols_A = i_cols;
119
120
        i_cols_Ap = i_cols_p;
        i_cols_As = i_cols_s;
121
        i_cols_Af = i_cols_f - ny;
122
123
124
125
126
127
128
129
130
        for i = 1:order+1
            i_w_p = 1;
            for j = 1:nnodes^(i-1)
                innovation = exo_simul;
                if i > 1
                    innovation(i+1,:) = nodes(mod(j-1,nnodes)+1,:);
                end
                if i <= order
                    for k=1:nnodes
131
132
                        y = [Y(i_cols_p,i_w_p);
                             Y(i_cols_s,j);
133
                             Y(i_cols_f,(j-1)*nnodes+k)];
134
135
136
137
138
                        [d1,jacobian] = dynamic_model(y,innovation,params,steady_state,i+1);
                        if i == 1
                            % in first period we don't keep track of
                            % predetermined variables
                            i_cols_A = [i_cols_As - ny; i_cols_Af];
139
                            A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_1);
140
141
                        else
                            i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af];
142
                            A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_j);
143
                        end
144
                        res(:,i,j) = res(:,i,j)+weights(k)*d1;
145
                        i_cols_Af = i_cols_Af + ny;
146
                    end
147
148
149
150
151
152
153
154
155
156
157
158
159
                else
                    y = [Y(i_cols_p,i_w_p);
                         Y(i_cols_s,j);
                         Y(i_cols_f,j)];
                    [d1,jacobian] = dynamic_model(y,innovation,params,steady_state,i+1);
                    if i == 1
                        % in first period we don't keep track of
                        % predetermined variables
                        i_cols_A = [i_cols_As - ny; i_cols_Af];
                        A1(i_rows,i_cols_A) = jacobian(:,i_cols_1);
                    else
                        i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af];
                        A1(i_rows,i_cols_A) = jacobian(:,i_cols_j);
160
                    end
161
162
                    res(:,i,j) = d1;
                    i_cols_Af = i_cols_Af + ny;
163
                end
164
165
166
                i_rows = i_rows + ny;
                if mod(j,nnodes) == 0
                    i_w_p = i_w_p + 1;
167
                end
168
169
170
171
172
                if i > 1
                    if mod(j,nnodes) == 0
                        i_cols_Ap = i_cols_Ap + ny;
                    end
                    i_cols_As = i_cols_As + ny;
173
                end
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
            end
            i_cols_p = i_cols_p + ny;
            i_cols_s = i_cols_s + ny;
            i_cols_f = i_cols_f + ny;
        end
        nzA = cell(periods,world_nbr);
        parfor j=1:world_nbr
            i_rows_y = find(lead_lag_incidence')+(order+1)*ny;
            offset_c = ny*(sum(nnodes.^(0:order-1),2)+j-1);
            offset_r = (j-1)*ny;
            for i=order+2:periods
                [d1,jacobian] = dynamic_model(Y(i_rows_y,j), ...
                                              exo_simul,params, ...
                                              steady_state,i+1);
                if i == periods
                    [ir,ic,v] = find(jacobian(:,i_cols_T));
                else
                    [ir,ic,v] = find(jacobian(:,i_cols_j));
Stéphane Adjemian's avatar
Stéphane Adjemian committed
192
                end
193
194
195
196
197
                nzA{i,j} = [offset_r+ir,offset_c+icA(ic), v]';
                res(:,i,j) = d1;
                i_rows_y = i_rows_y + ny;
                offset_c = offset_c + world_nbr*ny;
                offset_r = offset_r + world_nbr*ny;
198
199
            end
        end
200
        err = max(abs(res(i_upd_r)));
201
        if err < tolerance
202
            stop = 1;
203
            if verbose
204
205
206
207
208
209
210
                fprintf('\n') ;
                disp([' Total time of simulation        :' num2str(etime(clock,h1))]) ;
                fprintf('\n') ;
                disp([' Convergency obtained.']) ;
                fprintf('\n') ;
            end
            flag = 0;% Convergency obtained.
211
            endo_simul = reshape(Y(:,1),ny,periods+2);%Y(ny+(1:ny),1);
212
213
                                                      %            figure;plot(Y(16:ny:(periods+2)*ny,:))
                                                      %            pause
214
215
            break
        end
216
217
218
219
        A2 = [nzA{:}]';
        A = [A1; sparse(A2(:,1),A2(:,2),A2(:,3),ny*(periods-order-1)*world_nbr,dimension)];
        dy = -A\res(i_upd_r);
        Y(i_upd_y) =   Y(i_upd_y) + dy;
220
221
222
    end

    if ~stop
223
        if verbose
224
225
226
            fprintf('\n') ;
            disp(['     Total time of simulation        :' num2str(etime(clock,h1))]) ;
            fprintf('\n') ;
227
            disp(['WARNING : maximum number of iterations is reached (modify options_.simul.maxit).']) ;
228
229
230
231
232
            fprintf('\n') ;
        end
        flag = 1;% more iterations are needed.
        endo_simul = 1;
    end
233
    if verbose
234
235
        disp (['-----------------------------------------------------']) ;
    end