1function [dr, info] = dyn_first_order_solver(jacobia, DynareModel, dr, DynareOptions, task)
2
3% Computes the first order reduced form of a DSGE model.
4%
5% INPUTS
6% - jacobia       [double]    matrix, the jacobian of the dynamic model.
7% - DynareModel   [struct]    Matlab's structre describing the model, M_ global.
8% - dr            [struct]    Matlab's structure describing the reduced form model.
9% - DynareOptions [struct]    Matlab's structure containing the current state of the options, oo_ global.
10% - task          [integer]   scalar, if task = 0 then decision rules are computed and if task = 1 then only eigenvales are computed.
11%
12% OUTPUTS
13% - dr            [struct]    Matlab's structure describing the reduced form model.
14% - info          [integer]   scalar, error code. Possible values are:
15%
16%                                     info=0 -> no error,
17%                                     info=1 -> the model doesn't determine the current variables uniquely,
18%                                     info=2 -> mjdgges dll returned an error,
19%                                     info=3 -> Blanchard and Kahn conditions are not satisfied: no stable equilibrium,
20%                                     info=4 -> Blanchard and Kahn conditions are not satisfied: indeterminacy,
21%                                     info=5 -> Blanchard and Kahn conditions are not satisfied: indeterminacy due to rank failure,
22%                                     info=7 -> One of the eigenvalues is close to 0/0 (infinity of complex solutions)
23
24% Copyright (C) 2001-2018 Dynare Team
25%
26% This file is part of Dynare.
27%
28% Dynare is free software: you can redistribute it and/or modify
29% it under the terms of the GNU General Public License as published by
30% the Free Software Foundation, either version 3 of the License, or
31% (at your option) any later version.
32%
33% Dynare is distributed in the hope that it will be useful,
34% but WITHOUT ANY WARRANTY; without even the implied warranty of
35% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
36% GNU General Public License for more details.
37%
38% You should have received a copy of the GNU General Public License
39% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
40
41persistent reorder_jacobian_columns innovations_idx index_s index_m index_c
42persistent index_p row_indx index_0m index_0p k1 k2 state_var
43persistent ndynamic nstatic nfwrd npred nboth nd nsfwrd n_current index_d
44persistent index_e index_d1 index_d2 index_e1 index_e2 row_indx_de_1
45persistent row_indx_de_2 cols_b
46
47
48if ~nargin
49    if nargout
50        error('dyn_first_order_solver:: Initialization mode returns zero argument!')
51    end
52    reorder_jacobian_columns = [];
53    return
54end
55
56exo_nbr = DynareModel.exo_nbr;
57
58if isempty(reorder_jacobian_columns)
59
60    maximum_lag = DynareModel.maximum_endo_lag;
61    kstate   = dr.kstate;
62    nfwrd    = DynareModel.nfwrd;
63    nboth    = DynareModel.nboth;
64    npred    = DynareModel.npred;
65    nstatic  = DynareModel.nstatic;
66    ndynamic = DynareModel.ndynamic;
67    nsfwrd   = DynareModel.nsfwrd;
68    n        = DynareModel.endo_nbr;
69
70    k1 = 1:(npred+nboth);
71    k2 = 1:(nfwrd+nboth);
72
73    order_var = dr.order_var;
74    nd = size(kstate,1);
75    lead_lag_incidence = DynareModel.lead_lag_incidence;
76    nz = nnz(lead_lag_incidence);
77
78    lead_id = find(lead_lag_incidence(maximum_lag+2,:));
79    lead_idx = lead_lag_incidence(maximum_lag+2,lead_id);
80    if maximum_lag
81        lag_id = find(lead_lag_incidence(1,:));
82        lag_idx = lead_lag_incidence(1,lag_id);
83        static_id = find((lead_lag_incidence(1,:) == 0) & ...
84                         (lead_lag_incidence(3,:) == 0));
85    else
86        lag_id = [];
87        lag_idx = [];
88        static_id = find(lead_lag_incidence(2,:)==0);
89    end
90
91    both_id = intersect(lead_id,lag_id);
92    if maximum_lag
93        no_both_lag_id = setdiff(lag_id,both_id);
94    else
95        no_both_lag_id = lag_id;
96    end
97    innovations_idx = nz+(1:exo_nbr);
98    state_var  = [no_both_lag_id, both_id];
99
100    index_c  = nonzeros(lead_lag_incidence(maximum_lag+1,:));             % Index of all endogenous variables present at time=t
101    n_current = length(index_c);
102
103    index_s  = npred+nboth+(1:nstatic);     % Index of all static
104                                            % endogenous variables
105                                            % present at time=t
106    indexi_0 = npred+nboth;
107
108    npred0 = nnz(lead_lag_incidence(maximum_lag+1,no_both_lag_id));
109    index_0m = indexi_0+nstatic+(1:npred0);
110    nfwrd0 = nnz(lead_lag_incidence(maximum_lag+1,lead_id));
111    index_0p = indexi_0+nstatic+npred0+(1:nfwrd0);
112    index_m  = 1:(npred+nboth);
113    index_p  = npred+nboth+n_current+(1:nfwrd+nboth);
114    row_indx_de_1 = 1:ndynamic;
115    row_indx_de_2 = ndynamic+(1:nboth);
116    row_indx = nstatic+row_indx_de_1;
117    index_d = [index_0m index_p];
118    llx = lead_lag_incidence(:,order_var);
119    index_d1 = [find(llx(maximum_lag+1,nstatic+(1:npred))), npred+nboth+(1:nsfwrd) ];
120    index_d2 = npred+(1:nboth);
121
122    index_e = [index_m index_0p];
123    index_e1 = [1:npred+nboth, npred+nboth+find(llx(maximum_lag+1,nstatic+npred+(1: ...
124                                                      nsfwrd)))];
125    index_e2 = npred+nboth+(1:nboth);
126
127    [~,cols_b] = find(lead_lag_incidence(maximum_lag+1, order_var));
128
129    reorder_jacobian_columns = [nonzeros(lead_lag_incidence(:,order_var)'); ...
130                        nz+(1:exo_nbr)'];
131end
132
133dr.ghx = [];
134dr.ghu = [];
135
136dr.state_var = state_var;
137
138jacobia = jacobia(:,reorder_jacobian_columns);
139
140if nstatic > 0
141    [Q, ~] = qr(jacobia(:,index_s));
142    aa = Q'*jacobia;
143else
144    aa = jacobia;
145end
146
147A = aa(:,index_m);  % Jacobain matrix for lagged endogeneous variables
148B(:,cols_b) = aa(:,index_c);  % Jacobian matrix for contemporaneous endogeneous variables
149C = aa(:,index_p);  % Jacobain matrix for led endogeneous variables
150
151info = 0;
152if task ~= 1 && (DynareOptions.dr_cycle_reduction || DynareOptions.dr_logarithmic_reduction)
153    if n_current < DynareModel.endo_nbr
154        if DynareOptions.dr_cycle_reduction
155            error(['The cycle reduction algorithme can''t be used when the ' ...
156                   'coefficient matrix for current variables isn''t invertible'])
157        elseif DynareOptions.dr_logarithmic_reduction
158            error(['The logarithmic reduction algorithme can''t be used when the ' ...
159                   'coefficient matrix for current variables isn''t invertible'])
160        end
161    end
162    if DynareOptions.gpu
163        gpuArray(A1);
164        gpuArray(B1);
165        gpuArray(C1);
166    end
167    A1 = [aa(row_indx,index_m ) zeros(ndynamic,nfwrd)];
168    B1 = [aa(row_indx,index_0m) aa(row_indx,index_0p) ];
169    C1 = [zeros(ndynamic,npred) aa(row_indx,index_p)];
170    if DynareOptions.dr_cycle_reduction
171        [ghx, info] = cycle_reduction(A1, B1, C1, DynareOptions.dr_cycle_reduction_tol);
172    else
173        [ghx, info] = logarithmic_reduction(C1, B1, A1, DynareOptions.dr_logarithmic_reduction_tol, DynareOptions.dr_logarithmic_reduction_maxiter);
174    end
175    if info
176        % cycle_reduction or logarithmic redution failed and set info
177        return
178    end
179    ghx = ghx(:,index_m);
180    hx = ghx(1:npred+nboth,:);
181    gx = ghx(1+npred:end,:);
182else
183    D = zeros(ndynamic+nboth);
184    E = zeros(ndynamic+nboth);
185    D(row_indx_de_1,index_d1) = aa(row_indx,index_d);
186    D(row_indx_de_2,index_d2) = eye(nboth);
187    E(row_indx_de_1,index_e1) = -aa(row_indx,index_e);
188    E(row_indx_de_2,index_e2) = eye(nboth);
189
190    [err, ss, tt, w, sdim, dr.eigval, info1] = mjdgges(E, D, DynareOptions.qz_criterium, DynareOptions.qz_zero_threshold);
191    mexErrCheck('mjdgges', err);
192
193    if info1
194        if info1 == -30
195            % one eigenvalue is close to 0/0
196            info(1) = 7;
197        else
198            info(1) = 2;
199            info(2) = info1;
200            info(3) = size(E,2);
201        end
202        return
203    end
204
205    dr.sdim = sdim;                      % Number of stable eigenvalues.
206    dr.edim = length(dr.eigval)-sdim;    % Number of exposive eigenvalues.
207
208    nba = nd-sdim;
209
210    if task==1
211        if rcond(w(npred+nboth+1:end,npred+nboth+1:end)) < 1e-9
212            dr.full_rank = 0;
213        else
214            dr.full_rank = 1;
215        end
216    end
217
218    if nba ~= nsfwrd
219        temp = sort(abs(dr.eigval));
220        if nba > nsfwrd
221            temp = temp(nd-nba+1:nd-nsfwrd)-1-DynareOptions.qz_criterium;
222            info(1) = 3;
223        elseif nba < nsfwrd
224            temp = temp(nd-nsfwrd+1:nd-nba)-1-DynareOptions.qz_criterium;
225            info(1) = 4;
226        end
227        info(2) = temp'*temp;
228        return
229    end
230
231    if task==1, return, end
232
233    %First order approximation
234    indx_stable_root = 1: (nd - nsfwrd);         %=> index of stable roots
235    indx_explosive_root = npred + nboth + 1:nd;  %=> index of explosive roots
236                                                 % derivatives with respect to dynamic state variables
237                                                 % forward variables
238    Z = w';
239    Z11 = Z(indx_stable_root,    indx_stable_root);
240    Z21  = Z(indx_explosive_root, indx_stable_root);
241    Z22  = Z(indx_explosive_root, indx_explosive_root);
242    opts.TRANSA = false; % needed by Octave 4.0.0
243    [minus_gx,rc] = linsolve(Z22,Z21,opts);
244    if rc < 1e-9
245        % Z22 is near singular
246        info(1) = 5;
247        info(2) = -log(rc);
248        return
249    end
250    gx  = -minus_gx;
251    % predetermined variables
252    opts.UT = true;
253    opts.TRANSA = true;
254    hx1 = linsolve(tt(indx_stable_root, indx_stable_root),Z11,opts)';
255    opts.UT = false;      % needed by Octave 4.0.0
256    opts.TRANSA = false;  % needed by Octave 4.0.0
257    hx2 = linsolve(Z11,ss(indx_stable_root, indx_stable_root)',opts)';
258    hx =  hx1*hx2;
259    ghx = [hx(k1,:); gx(k2(nboth+1:end),:)];
260end
261
262if nstatic > 0
263    B_static = B(:,1:nstatic);  % submatrix containing the derivatives w.r. to static variables
264else
265    B_static = [];
266end
267%static variables, backward variable, mixed variables and forward variables
268B_pred = B(:,nstatic+1:nstatic+npred+nboth);
269B_fyd = B(:,nstatic+npred+nboth+1:end);
270
271% static variables
272if nstatic > 0
273    temp = - C(1:nstatic,:)*gx*hx;
274    b(:,cols_b) = aa(:,index_c);
275    b10 = b(1:nstatic, 1:nstatic);
276    b11 = b(1:nstatic, nstatic+1:end);
277    temp(:,index_m) = temp(:,index_m)-A(1:nstatic,:);
278    temp = b10\(temp-b11*ghx);
279    ghx = [temp; ghx];
280    temp = [];
281end
282
283A_ = real([B_static C*gx+B_pred B_fyd]); % The state_variable of the block are located at [B_pred B_both]
284
285if exo_nbr
286    if nstatic > 0
287        fu = Q' * jacobia(:,innovations_idx);
288    else
289        fu = jacobia(:,innovations_idx);
290    end
291
292    ghu = - A_ \ fu;
293else
294    ghu = [];
295end
296
297dr.ghx = ghx;
298dr.ghu = ghu;
299
300if DynareOptions.aim_solver ~= 1
301    % Necessary when using Sims' routines for QZ
302    dr.ghx = real(ghx);
303    dr.ghu = real(ghu);
304    hx = real(hx);
305end
306
307% non-predetermined variables
308dr.gx = gx;
309%predetermined (endogenous state) variables, square transition matrix
310dr.Gy = hx;
311