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