dyn_first_order_solver.m 10.6 KB
Newer Older
1
2
function [dr,info] = dyn_first_order_solver(jacobia,DynareModel,dr,DynareOptions,task)

3
%@info:
4
%! @deftypefn {Function File} {[@var{dr},@var{info}] =} dyn_first_order_solver (@var{jacobia},@var{DynareModel},@var{dr},@var{DynareOptions},@var{task})
5
6
7
8
9
10
11
12
13
%! @anchor{dyn_first_order_solver}
%! @sp 1
%! Computes the first order reduced form of the DSGE model
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item jacobia
%! Matrix containing the Jacobian of the model
14
%! @item DynareModel
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
%! Matlab's structure describing the model (initialized by @code{dynare}).
%! @item dr
%! Matlab's structure describing the reduced form solution of the model.
%! @item qz_criterium
%! Double containing the criterium to separate explosive from stable eigenvalues
%! @end table
%! @sp 2
%! @strong{Outputs}
%! @sp 1
%! @table @ @var
%! @item dr
%! Matlab's structure describing the reduced form solution of the model.
%! @item info
%! Integer scalar, error code.
%! @sp 1
%! @table @ @code
%! @item info==0
%! No error.
%! @item info==1
%! The model doesn't determine the current variables uniquely.
%! @item info==2
%! MJDGGES returned an error code.
%! @item info==3
%! Blanchard & Kahn conditions are not satisfied: no stable equilibrium.
%! @item info==4
%! Blanchard & Kahn conditions are not satisfied: indeterminacy.
%! @item info==5
%! Blanchard & Kahn conditions are not satisfied: indeterminacy due to rank failure.
%! @item info==7
44
%! One of the generalized eigenvalues is close to 0/0
45
%! @end table
MichelJuillard's avatar
MichelJuillard committed
46
%! @end table
47
48
49
%! @end deftypefn
%@eod:

Sébastien Villemot's avatar
Sébastien Villemot committed
50
% Copyright (C) 2001-2012 Dynare Team
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
%
% 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/>.
66

67
68
persistent reorder_jacobian_columns innovations_idx index_s index_m index_c
persistent index_p row_indx index_0m index_0p k1 k2 state_var
69
persistent ndynamic nstatic nfwrd npred nboth nd nsfwrd n_current index_d 
70
71
persistent index_e index_d1 index_d2 index_e1 index_e2 row_indx_de_1 
persistent row_indx_de_2 cols_b
72
73
74


if ~nargin
75
76
77
    if nargout
        error('dyn_first_order_solver:: Initialization mode returns zero argument!')
    end
78
79
80
81
    reorder_jacobian_columns = [];
    return
end

82
exo_nbr = DynareModel.exo_nbr;
83

84
85
86
if isempty(reorder_jacobian_columns)
    
    maximum_lag = DynareModel.maximum_endo_lag;
87
    kstate   = dr.kstate;
88
89
90
91
92
93
    nfwrd    = DynareModel.nfwrd;
    nboth    = DynareModel.nboth;
    npred    = DynareModel.npred;
    nstatic  = DynareModel.nstatic;
    ndynamic = DynareModel.ndynamic;
    nsfwrd   = DynareModel.nsfwrd;
94
    n        = DynareModel.endo_nbr;
95
96
97
98

    k1 = 1:(npred+nboth);
    k2 = 1:(nfwrd+nboth);

99
100
    order_var = dr.order_var;
    nd = size(kstate,1);
101
    lead_lag_incidence = DynareModel.lead_lag_incidence;
102
103
    nz = nnz(lead_lag_incidence);

104
105
106
107
108
109
110
111
112
113
114
    lead_id = find(lead_lag_incidence(maximum_lag+2,:));
    lead_idx = lead_lag_incidence(maximum_lag+2,lead_id);
    if maximum_lag
        lag_id = find(lead_lag_incidence(1,:));
        lag_idx = lead_lag_incidence(1,lag_id);
        static_id = find((lead_lag_incidence(1,:) == 0) & ...
                         (lead_lag_incidence(3,:) == 0));
    else
        lag_id = [];
        lag_idx = [];
        static_id = find(lead_lag_incidence(2,:)==0);
115
116
    end

117
118
119
120
    both_id = intersect(lead_id,lag_id);
    if maximum_lag
        no_both_lag_id = setdiff(lag_id,both_id);
    else
121
        no_both_lag_id = lag_id;
122
123
    end
    innovations_idx = nz+(1:exo_nbr);
124
    state_var  = [no_both_lag_id, both_id];
125

126
127
    index_c  = nonzeros(lead_lag_incidence(maximum_lag+1,:));             % Index of all endogenous variables present at time=t
    n_current = length(index_c);
128

129
130
131
132
    index_s  = npred+nboth+(1:nstatic);     % Index of all static
                                            % endogenous variables
                                            % present at time=t
    indexi_0 = npred+nboth;
133

134
135
136
137
138
    npred0 = nnz(lead_lag_incidence(maximum_lag+1,no_both_lag_id));
    index_0m = indexi_0+nstatic+(1:npred0);
    nfwrd0 = nnz(lead_lag_incidence(2,lead_id));
    index_0p = indexi_0+nstatic+npred0+(1:nfwrd0);
    index_m  = 1:(npred+nboth);
139
    index_p  = npred+nboth+n_current+(1:nfwrd+nboth);
140
    row_indx_de_1 = 1:ndynamic;
141
    row_indx_de_2 = ndynamic+(1:nboth);
142
143
144
    row_indx = nstatic+row_indx_de_1;
    index_d = [index_0m index_p];
    llx = lead_lag_incidence(:,order_var);
145
    index_d1 = [find(llx(maximum_lag+1,nstatic+(1:npred))), npred+nboth+(1:nsfwrd) ];
146
    index_d2 = npred+(1:nboth);
147
148

    index_e = [index_m index_0p];
149
    index_e1 = [1:npred+nboth, npred+nboth+find(llx(maximum_lag+1,nstatic+npred+(1: ...
150
                                                      nsfwrd)))];
151
    index_e2 = npred+nboth+(1:nboth);
152
153
154
155
156
157
    
    [junk,cols_b] = find(lead_lag_incidence(maximum_lag+1, order_var));

    reorder_jacobian_columns = [nonzeros(lead_lag_incidence(:,order_var)'); ...
                        nz+(1:exo_nbr)'];
end
158
159
160
161

dr.ghx = [];
dr.ghu = [];

162
163
dr.state_var = state_var;

164
165
166
167
168
169
170
171
172
173
jacobia = jacobia(:,reorder_jacobian_columns);

if nstatic > 0
    [Q, junk] = qr(jacobia(:,index_s));
    aa = Q'*jacobia;
else
    aa = jacobia;
end

A = aa(:,index_m);  % Jacobain matrix for lagged endogeneous variables
174
B(:,cols_b) = aa(:,index_c);  % Jacobian matrix for contemporaneous endogeneous variables
175
176
C = aa(:,index_p);  % Jacobain matrix for led endogeneous variables

177
info = 0;
178
if task ~= 1 && (DynareOptions.dr_cycle_reduction || DynareOptions.dr_logarithmic_reduction)
179
180
181
    if n_current < DynareModel.endo_nbr
        if DynareOptions.dr_cycle_reduction
            error(['The cycle reduction algorithme can''t be used when the ' ...
MichelJuillard's avatar
MichelJuillard committed
182
               'coefficient matrix for current variables isn''t invertible'])
183
184
        elseif DynareOptions.dr_logarithmic_reduction
            error(['The logarithmic reduction algorithme can''t be used when the ' ...
MichelJuillard's avatar
MichelJuillard committed
185
                   'coefficient matrix for current variables isn''t invertible'])
186
        end
187
188
189
190
191
192
    end  
    if DynareOptions.gpu
        gpuArray(A1);
        gpuArray(B1);
        gpuArray(C1);
    end
193
194
195
    A1 = [aa(row_indx,index_m ) zeros(ndynamic,nfwrd)];
    B1 = [aa(row_indx,index_0m) aa(row_indx,index_0p) ];
    C1 = [zeros(ndynamic,npred) aa(row_indx,index_p)];
196
    if DynareOptions.dr_cycle_reduction == 1
MichelJuillard's avatar
MichelJuillard committed
197
        [ghx, info] = cycle_reduction(A1, B1, C1, DynareOptions.dr_cycle_reduction_tol);
198
    else
MichelJuillard's avatar
MichelJuillard committed
199
200
201
202
203
        [ghx, info] = logarithmic_reduction(C1, B1, A1, DynareOptions.dr_logarithmic_reduction_tol, DynareOptions.dr_logarithmic_reduction_maxiter);
    end
    if info
        % cycle_reduction or logarithmic redution failed and set info
        return
204
    end
205
206
207
    ghx = ghx(:,index_m);
    hx = ghx(1:npred+nboth,:);
    gx = ghx(1+npred:end,:);
MichelJuillard's avatar
MichelJuillard committed
208
else
209
210
211
212
213
214
215
    D = zeros(ndynamic+nboth);
    E = zeros(ndynamic+nboth);
    D(row_indx_de_1,index_d1) = aa(row_indx,index_d);
    D(row_indx_de_2,index_d2) = eye(nboth);
    E(row_indx_de_1,index_e1) = -aa(row_indx,index_e);
    E(row_indx_de_2,index_e2) = eye(nboth);
    
216
    [err, ss, tt, w, sdim, dr.eigval, info1] = mjdgges(E, D, DynareOptions.qz_criterium, DynareOptions.qz_zero_threshold);
217
218
    mexErrCheck('mjdgges', err);

MichelJuillard's avatar
MichelJuillard committed
219
220
    if info1
        if info1 == -30
221
            % one eigenvalue is close to 0/0
222
223
224
225
            info(1) = 7;
        else
            info(1) = 2;
            info(2) = info1;
226
            info(3) = size(E,2);
227
228
229
230
231
232
233
        end
        return
    end

    nba = nd-sdim;

    if task == 1
234
235
236
237
238
        if rcond(w(npred+nboth+1:end,npred+nboth+1:end)) < 1e-9
            dr.full_rank = 0;
        else
            dr.full_rank = 1;
        end            
239
240
241
        return
    end

242
    if nba ~= nsfwrd
243
        temp = sort(abs(dr.eigval));
244
245
        if nba > nsfwrd
            temp = temp(nd-nba+1:nd-nsfwrd)-1-DynareOptions.qz_criterium;
246
            info(1) = 3;
247
248
        elseif nba < nsfwrd;
            temp = temp(nd-nsfwrd+1:nd-nba)-1-DynareOptions.qz_criterium;
249
250
251
252
253
254
            info(1) = 4;
        end
        info(2) = temp'*temp;
        return
    end

255
    %First order approximation
256
    indx_stable_root = 1: (nd - nsfwrd);     %=> index of stable roots
257
258
259
260
    indx_explosive_root = npred + nboth + 1:nd;  %=> index of explosive roots
                                                 % derivatives with respect to dynamic state variables
                                                 % forward variables
    Z = w';
MichelJuillard's avatar
MichelJuillard committed
261
    Z11 = Z(indx_stable_root,    indx_stable_root);
262
263
    Z21  = Z(indx_explosive_root, indx_stable_root);
    Z22  = Z(indx_explosive_root, indx_explosive_root);
264
265
    opts.TRANSA = false; % needed by Octave 4.0.0
    [minus_gx,rc] = linsolve(Z22,Z21,opts);
MichelJuillard's avatar
MichelJuillard committed
266
267
    if rc < 1e-9
        % Z22 is near singular
268
        info(1) = 5;
MichelJuillard's avatar
MichelJuillard committed
269
        info(2) = -log(rc);
270
271
        return;
    end
MichelJuillard's avatar
MichelJuillard committed
272
    gx  = -minus_gx;
273
    % predetermined variables
MichelJuillard's avatar
MichelJuillard committed
274
275
276
    opts.UT = true;
    opts.TRANSA = true;
    hx1 = linsolve(tt(indx_stable_root, indx_stable_root),Z11,opts)';
277
278
    opts.UT = false;      % needed by Octave 4.0.0
    opts.TRANSA = false;  % needed by Octave 4.0.0
279
    hx2 = linsolve(Z11,ss(indx_stable_root, indx_stable_root)',opts)';
MichelJuillard's avatar
MichelJuillard committed
280
    hx =  hx1*hx2;
281
282
    ghx = [hx(k1,:); gx(k2(nboth+1:end),:)];
end
283

284
285
286
287
288
289
290
291
if nstatic > 0
    B_static = B(:,1:nstatic);  % submatrix containing the derivatives w.r. to static variables
else
    B_static = [];
end;
%static variables, backward variable, mixed variables and forward variables
B_pred = B(:,nstatic+1:nstatic+npred+nboth);
B_fyd = B(:,nstatic+npred+nboth+1:end);
292

293
294
295
296
297
298
299
300
301
302
303
% static variables
if nstatic > 0
    temp = - C(1:nstatic,:)*gx*hx;
    b(:,cols_b) = aa(:,index_c);
    b10 = b(1:nstatic, 1:nstatic);
    b11 = b(1:nstatic, nstatic+1:end);
    temp(:,index_m) = temp(:,index_m)-A(1:nstatic,:);
    temp = b10\(temp-b11*ghx);
    ghx = [temp; ghx];
    temp = [];
end
304

305
A_ = real([B_static C*gx+B_pred B_fyd]); % The state_variable of the block are located at [B_pred B_both]
306

307
308
309
if exo_nbr
    if nstatic > 0
        fu = Q' * jacobia(:,innovations_idx);
310
    else
311
        fu = jacobia(:,innovations_idx);
312
    end;
313
314
315
316
317
318

    ghu = - A_ \ fu;
else
    ghu = [];
end;

319
320
321
322
323
324
325
326
327
dr.ghx = ghx;
dr.ghu = ghu;

if DynareOptions.aim_solver ~= 1 && DynareOptions.use_qzdiv
    % Necessary when using Sims' routines for QZ
    dr.ghx = real(ghx);
    dr.ghu = real(ghu);
    hx = real(hx);
end
328

329
330
331
% non-predetermined variables
dr.gx = gx;
%predetermined (endogenous state) variables, square transition matrix
332
dr.Gy = hx;