1
2
3
4
5
6
7
8
9
10
11
12
13
14
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
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
|
function dr = dyn_second_order_solver(jacobia,hessian_mat,dr,M_,threads_BC)
%@info:
%! @deftypefn {Function File} {@var{dr} =} dyn_second_order_solver (@var{jacobia},@var{hessian_mat},@var{dr},@var{M_},@var{threads_BC})
%! @anchor{dyn_second_order_solver}
%! @sp 1
%! Computes the second order reduced form of the DSGE model, for details please refer to
%! * Juillard and Kamenik (2004): Solving Stochastic Dynamic Equilibrium Models: A k-Order Perturbation Approach
%! * Kamenik (2005) - Solving SDGE Models: A New Algorithm for the Sylvester Equation
%! Note that this function makes use of the fact that Dynare internally transforms the model
%! so that there is only one lead and one lag on endogenous variables and, in the case of a stochastic model,
%! no leads/lags on exogenous variables. See the manual for more details.
% Auxiliary variables
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item jacobia
%! Matrix containing the Jacobian of the model
%! @item hessian_mat
%! Matrix containing the second order derivatives of the model
%! @item dr
%! Matlab's structure describing the reduced form solution of the model.
%! @item M_
%! Matlab's structure describing the model (initialized by @code{dynare}).
%! @item threads_BC
%! Integer controlling number of threads in sparse_hessian_times_B_kronecker_C
%! @end table
%! @sp 2
%! @strong{Outputs}
%! @sp 1
%! @table @ @var
%! @item dr
%! Matlab's structure describing the reduced form solution of the model.
%! @end table
%! @end deftypefn
%@eod:
% Copyright © 2001-2020 Dynare Team
%
% 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 <https://www.gnu.org/licenses/>.
dr.ghxx = [];
dr.ghuu = [];
dr.ghxu = [];
dr.ghs2 = [];
k1 = nonzeros(M_.lead_lag_incidence(:,dr.order_var)');
kk1 = [k1; length(k1)+(1:M_.exo_nbr+M_.exo_det_nbr)'];
nk = size(kk1,1);
kk2 = reshape(1:nk^2,nk,nk);
ic = [ M_.nstatic+(1:M_.nspred) ]';
klag = M_.lead_lag_incidence(1,dr.order_var); %columns are in DR order
kcurr = M_.lead_lag_incidence(2,dr.order_var); %columns are in DR order
klead = M_.lead_lag_incidence(3,dr.order_var); %columns are in DR order
%% ghxx
A = zeros(M_.endo_nbr,M_.endo_nbr);
A(:,kcurr~=0) = jacobia(:,nonzeros(kcurr));
A(:,ic) = A(:,ic) + jacobia(:,nonzeros(klead))*dr.ghx(klead~=0,:);
B = zeros(M_.endo_nbr,M_.endo_nbr);
B(:,M_.nstatic+M_.npred+1:end) = jacobia(:,nonzeros(klead));
C = dr.ghx(ic,:);
zx = [eye(length(ic));
dr.ghx(kcurr~=0,:);
dr.ghx(klead~=0,:)*dr.ghx(ic,:);
zeros(M_.exo_nbr,length(ic));
zeros(M_.exo_det_nbr,length(ic))];
zu = [zeros(length(ic),M_.exo_nbr);
dr.ghu(kcurr~=0,:);
dr.ghx(klead~=0,:)*dr.ghu(ic,:);
eye(M_.exo_nbr);
zeros(M_.exo_det_nbr,M_.exo_nbr)];
rhs = sparse_hessian_times_B_kronecker_C(hessian_mat(:,kk2(kk1,kk1)),zx,threads_BC); %hessian_mat: reordering to DR order
rhs = -rhs;
dr.ghxx = gensylv(2,A,B,C,rhs);
%% ghxu
%rhs
rhs = sparse_hessian_times_B_kronecker_C(hessian_mat(:,kk2(kk1,kk1)),zx,zu,threads_BC); %hessian_mat: reordering to DR order
abcOut = A_times_B_kronecker_C(dr.ghxx, dr.ghx(ic,:), dr.ghu(ic,:));
rhs = -rhs-B*abcOut;
%lhs
dr.ghxu = A\rhs;
%% ghuu
%rhs
rhs = sparse_hessian_times_B_kronecker_C(hessian_mat(:,kk2(kk1,kk1)),zu,threads_BC); %hessian_mat: reordering to DR order
B1 = A_times_B_kronecker_C(B*dr.ghxx,dr.ghu(ic,:));
rhs = -rhs-B1;
%lhs
dr.ghuu = A\rhs;
%% ghs2
% derivatives of F with respect to forward variables
O1 = zeros(M_.endo_nbr,M_.nstatic);
O2 = zeros(M_.endo_nbr,M_.nfwrd);
LHS = zeros(M_.endo_nbr,M_.endo_nbr);
LHS(:,kcurr~=0) = jacobia(:,nonzeros(kcurr));
RHS = zeros(M_.endo_nbr,M_.exo_nbr^2);
E = eye(M_.endo_nbr);
B1 = 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
RHS = RHS + jacobia(:,nonzeros(klead))*dr.ghuu(klead~=0,:)+B1;
% LHS
LHS = LHS + jacobia(:,nonzeros(klead))*(E(klead~=0,:)+[O1(klead~=0,:) dr.ghx(klead~=0,:) O2(klead~=0,:)]);
RHS = RHS*M_.Sigma_e(:);
dr.fuu = RHS;
RHS = -RHS;
dr.ghs2 = LHS\RHS;
|
