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