1function dr = dyn_second_order_solver(jacobia,hessian_mat,dr,M,threads_BC) 2 3%@info: 4%! @deftypefn {Function File} {@var{dr} =} dyn_second_order_solver (@var{jacobia},@var{hessian_mat},@var{dr},@var{M_},@var{threads_BC}) 5%! @anchor{dyn_second_order_solver} 6%! @sp 1 7%! Computes the second order reduced form of the DSGE model, for details please refer to 8%! * Juillard and Kamenik (2004): Solving Stochastic Dynamic Equilibrium Models: A k-Order Perturbation Approach 9%! * Kamenik (2005) - Solving SDGE Models: A New Algorithm for the Sylvester Equation 10%! Note that this function makes use of the fact that Dynare internally transforms the model 11%! so that there is only one lead and one lag on endogenous variables and, in the case of a stochastic model, 12%! no leads/lags on exogenous variables. See the manual for more details. 13% Auxiliary variables 14%! @sp 2 15%! @strong{Inputs} 16%! @sp 1 17%! @table @ @var 18%! @item jacobia 19%! Matrix containing the Jacobian of the model 20%! @item hessian_mat 21%! Matrix containing the second order derivatives of the model 22%! @item dr 23%! Matlab's structure describing the reduced form solution of the model. 24%! @item M_ 25%! Matlab's structure describing the model (initialized by @code{dynare}). 26%! @item threads_BC 27%! Integer controlling number of threads in sparse_hessian_times_B_kronecker_C 28%! @end table 29%! @sp 2 30%! @strong{Outputs} 31%! @sp 1 32%! @table @ @var 33%! @item dr 34%! Matlab's structure describing the reduced form solution of the model. 35%! @end table 36%! @end deftypefn 37%@eod: 38 39% Copyright (C) 2001-2019 Dynare Team 40% 41% This file is part of Dynare. 42% 43% Dynare is free software: you can redistribute it and/or modify 44% it under the terms of the GNU General Public License as published by 45% the Free Software Foundation, either version 3 of the License, or 46% (at your option) any later version. 47% 48% Dynare is distributed in the hope that it will be useful, 49% but WITHOUT ANY WARRANTY; without even the implied warranty of 50% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 51% GNU General Public License for more details. 52% 53% You should have received a copy of the GNU General Public License 54% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 55 56dr.ghxx = []; 57dr.ghuu = []; 58dr.ghxu = []; 59dr.ghs2 = []; 60 61k1 = nonzeros(M.lead_lag_incidence(:,dr.order_var)'); 62kk1 = [k1; length(k1)+(1:M.exo_nbr+M.exo_det_nbr)']; 63nk = size(kk1,1); 64kk2 = reshape(1:nk^2,nk,nk); 65ic = [ M.nstatic+(1:M.nspred) M.endo_nbr+(1:size(dr.ghx,2)-M.nspred) ]'; 66 67klag = M.lead_lag_incidence(1,dr.order_var); %columns are in DR order 68kcurr = M.lead_lag_incidence(2,dr.order_var); %columns are in DR order 69klead = M.lead_lag_incidence(3,dr.order_var); %columns are in DR order 70 71%% ghxx 72A = zeros(M.endo_nbr,M.endo_nbr); 73A(:,kcurr~=0) = jacobia(:,nonzeros(kcurr)); 74A(:,ic) = A(:,ic) + jacobia(:,nonzeros(klead))*dr.ghx(klead~=0,:); 75B = zeros(M.endo_nbr,M.endo_nbr); 76B(:,M.nstatic+M.npred+1:end) = jacobia(:,nonzeros(klead)); 77C = dr.ghx(ic,:); 78zx = [eye(length(ic)); 79 dr.ghx(kcurr~=0,:); 80 dr.ghx(klead~=0,:)*dr.ghx(ic,:); 81 zeros(M.exo_nbr,length(ic)); 82 zeros(M.exo_det_nbr,length(ic))]; 83zu = [zeros(length(ic),M.exo_nbr); 84 dr.ghu(kcurr~=0,:); 85 dr.ghx(klead~=0,:)*dr.ghu(ic,:); 86 eye(M.exo_nbr); 87 zeros(M.exo_det_nbr,M.exo_nbr)]; 88[rhs, err] = sparse_hessian_times_B_kronecker_C(hessian_mat(:,kk2(kk1,kk1)),zx,threads_BC); %hessian_mat: reordering to DR order 89mexErrCheck('sparse_hessian_times_B_kronecker_C', err); 90rhs = -rhs; 91[err, dr.ghxx] = gensylv(2,A,B,C,rhs); 92mexErrCheck('gensylv', err); 93 94 95%% ghxu 96%rhs 97[rhs, err] = sparse_hessian_times_B_kronecker_C(hessian_mat(:,kk2(kk1,kk1)),zx,zu,threads_BC); %hessian_mat: reordering to DR order 98mexErrCheck('sparse_hessian_times_B_kronecker_C', err); 99[abcOut,err] = A_times_B_kronecker_C(dr.ghxx, dr.ghx(ic,:), dr.ghu(ic,:)); 100mexErrCheck('A_times_B_kronecker_C', err); 101rhs = -rhs-B*abcOut; 102%lhs 103dr.ghxu = A\rhs; 104 105%% ghuu 106%rhs 107[rhs, err] = sparse_hessian_times_B_kronecker_C(hessian_mat(:,kk2(kk1,kk1)),zu,threads_BC); %hessian_mat: reordering to DR order 108mexErrCheck('sparse_hessian_times_B_kronecker_C', err); 109[B1, err] = A_times_B_kronecker_C(B*dr.ghxx,dr.ghu(ic,:)); 110mexErrCheck('A_times_B_kronecker_C', err); 111rhs = -rhs-B1; 112%lhs 113dr.ghuu = A\rhs; 114 115%% ghs2 116% derivatives of F with respect to forward variables 117O1 = zeros(M.endo_nbr,M.nstatic); 118O2 = zeros(M.endo_nbr,M.nfwrd); 119LHS = zeros(M.endo_nbr,M.endo_nbr); 120LHS(:,kcurr~=0) = jacobia(:,nonzeros(kcurr)); 121RHS = zeros(M.endo_nbr,M.exo_nbr^2); 122E = eye(M.endo_nbr); 123[B1, err] = sparse_hessian_times_B_kronecker_C(hessian_mat(:,kk2(nonzeros(klead),nonzeros(klead))), dr.ghu(klead~=0,:),threads_BC); %hessian_mat:focus only on forward variables and reorder to DR order 124mexErrCheck('sparse_hessian_times_B_kronecker_C', err); 125RHS = RHS + jacobia(:,nonzeros(klead))*dr.ghuu(klead~=0,:)+B1; 126% LHS 127LHS = LHS + jacobia(:,nonzeros(klead))*(E(klead~=0,:)+[O1(klead~=0,:) dr.ghx(klead~=0,:) O2(klead~=0,:)]); 128RHS = RHS*M.Sigma_e(:); 129dr.fuu = RHS; 130RHS = -RHS; 131dr.ghs2 = LHS\RHS; 132