dyn_second_order_solver.m 4.92 KB
Newer Older
1
function 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})
Sébastien Villemot's avatar
Sébastien Villemot committed
5
%! @anchor{dyn_second_order_solver}
6
%! @sp 1
7 8 9 10 11 12 13
%! 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
14 15 16 17 18 19
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item jacobia
%! Matrix containing the Jacobian of the model
20
%! @item hessian_mat
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
%! 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:
38

39
% Copyright (C) 2001-2019 Dynare Team
40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
%
% 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 <http://www.gnu.org/licenses/>.
55 56 57 58 59 60

dr.ghxx = [];
dr.ghuu = [];
dr.ghxu = [];
dr.ghs2 = [];

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
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) M.endo_nbr+(1:size(dr.ghx,2)-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, err] = sparse_hessian_times_B_kronecker_C(hessian_mat(:,kk2(kk1,kk1)),zx,threads_BC); %hessian_mat: reordering to DR order
89 90
mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
rhs = -rhs;
91
[err, dr.ghxx] = gensylv(2,A,B,C,rhs);
92 93
mexErrCheck('gensylv', err);

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
98
mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
99
[abcOut,err] = A_times_B_kronecker_C(dr.ghxx, dr.ghx(ic,:), dr.ghu(ic,:));
100
mexErrCheck('A_times_B_kronecker_C', err);
101
rhs = -rhs-B*abcOut;
102 103 104
%lhs
dr.ghxu = A\rhs;

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
108
mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
109
[B1, err] = A_times_B_kronecker_C(B*dr.ghxx,dr.ghu(ic,:));
110
mexErrCheck('A_times_B_kronecker_C', err);
111
rhs = -rhs-B1;
112 113 114
%lhs
dr.ghuu = A\rhs;

115
%% ghs2
116
% derivatives of F with respect to forward variables
117 118 119 120 121 122 123
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, 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
124
mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
125
RHS = RHS + jacobia(:,nonzeros(klead))*dr.ghuu(klead~=0,:)+B1;
126
% LHS
127 128
LHS = LHS + jacobia(:,nonzeros(klead))*(E(klead~=0,:)+[O1(klead~=0,:) dr.ghx(klead~=0,:) O2(klead~=0,:)]);
RHS = RHS*M.Sigma_e(:);
129 130 131
dr.fuu = RHS;
RHS = -RHS;
dr.ghs2 = LHS\RHS;