From a36287fc3dda2a7f02429e4f786b3adf358ddd7e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Fri, 29 Mar 2024 18:21:53 +0100
Subject: [PATCH] Convert stochastic solvers to sparse Jacobian representation
Ref. #1859
---
matlab/AIM/dynAIMsolver1.m | 57 ++-----
matlab/display_problematic_vars_Jacobian.m | 58 ++++---
matlab/model_diagnostics.m | 54 +++---
.../dyn_first_order_solver.m | 15 +-
.../dyn_second_order_solver.m | 47 +++---
matlab/stochastic_solver/stochastic_solvers.m | 156 +++++++++---------
6 files changed, 179 insertions(+), 208 deletions(-)
diff --git a/matlab/AIM/dynAIMsolver1.m b/matlab/AIM/dynAIMsolver1.m
index 78d770cd37..7c02fbf531 100644
--- a/matlab/AIM/dynAIMsolver1.m
+++ b/matlab/AIM/dynAIMsolver1.m
@@ -1,6 +1,6 @@
-function [dr,aimcode,rts]=dynAIMsolver1(jacobia_,M_,dr)
-% function [dr,aimcode]=dynAIMsolver1(jacobia_,M_,dr)
-% Maps Dynare jacobia to AIM 1st order model solver designed and developed by Gary ANderson
+function [dr, aimcode, rts] = dynAIMsolver1(g1, M_, dr)
+% function [dr, aimcode, rts] = dynAIMsolver1(g1, M_, dr)
+% Maps Dynare jacobian to AIM 1st order model solver designed and developed by Gary ANderson
% and derives the solution for gy=dr.hgx and gu=dr.hgu from the AIM outputs
% AIM System is given as a sum:
% i.e. for i=-$...+& SUM(Hi*xt+i)= £*zt, t = 0, . . . ,?
@@ -11,7 +11,7 @@ function [dr,aimcode,rts]=dynAIMsolver1(jacobia_,M_,dr)
% where [fy'-$... fy'i ... fy'+&]=[H-$... Hi ... H+&] and fu'= £
%
% INPUTS
-% jacobia_ [matrix] 1st order derivative of the model
+% g1 [matrix] sparse Jacobian of the dynamic model
% dr [matlab structure] Decision rules for stochastic simulations.
% M_ [matlab structure] Definition of the model.
%
@@ -46,7 +46,7 @@ function [dr,aimcode,rts]=dynAIMsolver1(jacobia_,M_,dr)
%
% GP July 2008
-% Copyright © 2008-2017 Dynare Team
+% Copyright © 2008-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -64,69 +64,36 @@ function [dr,aimcode,rts]=dynAIMsolver1(jacobia_,M_,dr)
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
aimcode=-1;
-neq= size(jacobia_,1); % no of equations
-lags=M_.maximum_endo_lag; % no of lags and leads
-leads=M_.maximum_endo_lead;
-klen=(leads+lags+1); % total lenght
-theAIM_H=zeros(neq, neq*klen); % alloc space
-lli=M_.lead_lag_incidence';
-% "sparse" the compact jacobia into AIM H aray of matrices
-% without exogenous shoc
-theAIM_H(:,find(lli(:)))=jacobia_(:,nonzeros(lli(:)));
+neq= size(g1, 1); % no of equations
condn = 1.e-10;%SPAmalg uses this in zero tests
uprbnd = 1 + 1.e-6;%allow unit roots
- % forward only models - AIM must have at least 1 lead and 1 lag.
-if lags ==0
- theAIM_H =[zeros(neq) theAIM_H];
- lags=1;
- klen=(leads+lags+1);
-end
-% backward looking only models
-if leads ==0
- theAIM_H =[theAIM_H zeros(neq)];
- leads=1;
- klen=(leads+lags+1);
-end
+
%disp('gensysToAMA:running ama');
try % try to run AIM
- [bb,rts,~,~,~,~,aimcode] =...
- SPAmalg(theAIM_H,neq, lags,leads,condn,uprbnd);
+ [bb, rts, ~, ~, ~, ~, aimcode] = SPAmalg(g1(:, 1:3*M_.endo_nbr), neq, 1, 1, condn, uprbnd);
catch
err = lasterror;
disp(['Dynare AIM Solver error:' sprintf('%s; ID:%s',err.message, err.identifier)]);
rethrow(lasterror);
end
if aimcode==1 %if OK
- col_order=[];
- for i =1:lags
- col_order=[((i-1)*neq)+dr.order_var' col_order];
- end
- bb_ord= bb(dr.order_var,col_order); % derive ordered gy
-
- % variables are present in the state space at the lag at which they
- % appear and at all smaller lags. The are ordered from smaller to
- % higher lag (reversed order of M_.lead_lag_incidence rows for lagged
- % variables)
- i_lagged_vars = flipud(cumsum(M_.lead_lag_incidence(1:M_.maximum_lag,dr.order_var),1))';
+ dr.ghx = bb(dr.order_var, dr.order_var(M_.nstatic+(1:M_.nspred)));
- dr.ghx= bb_ord(:,find(i_lagged_vars(:))); % get columns reported in
- % Dynare solution
if M_.exo_nbr % if there are exogenous shocks then derive gu for the shocks:
% get H0 and H+1=HM
% theH0= theAIM_H (:,M_.maximum_endo_lag*neq+1: (M_.maximum_endo_lag+1)*neq);
%theH0= theAIM_H (:,lags*neq+1: (lags+1)*neq);
% theHP= theAIM_H (:,(M_.maximum_endo_lag+1)*neq+1: (M_.maximum_endo_lag+2)*neq);
%theHP= theAIM_H (:,(lags+1)*neq+1: (lags+2)*neq);
- theAIM_Psi= - jacobia_(:, size(nonzeros(lli(:)))+1:end);%
%? = inv(H0 + H1B1)
%phi= (theH0+theHP*sparse(bb(:,(lags-1)*neq+1:end)))\eye( neq);
%AIM_ghu=phi*theAIM_Psi;
%dr.ghu =AIM_ghu(dr.order_var,:); % order gu
% Using AIM SPObstruct
- scof = SPObstruct(theAIM_H,bb,neq,lags,leads);
- scof1= scof(:,(lags)*neq+1:end);
+ scof = SPObstruct(g1(:, 1:3*M_.endo_nbr), bb, neq, 1, 1);
+ scof1 = scof(:, neq+1:end);
scof1= scof1(:,dr.order_var);
- dr.ghu =scof1\theAIM_Psi;
+ dr.ghu = -scof1 \ g1(:, 3*M_.endo_nbr+1:end);
else
dr.ghu = [];
end
diff --git a/matlab/display_problematic_vars_Jacobian.m b/matlab/display_problematic_vars_Jacobian.m
index f204b421a8..5105a3b301 100644
--- a/matlab/display_problematic_vars_Jacobian.m
+++ b/matlab/display_problematic_vars_Jacobian.m
@@ -1,5 +1,5 @@
-function []=display_problematic_vars_Jacobian(problemrow,problemcol,M_,x,type,caller_string)
-% []=display_problematic_vars_Jacobian(problemrow,problemcol,M_,x,caller_string)
+function display_problematic_vars_Jacobian(problemrow, problemcol, M_, x, type, caller_string)
+% display_problematic_vars_Jacobian(problemrow,problemcol,M_,x,caller_string)
% print the equation numbers and variables associated with problematic entries
% of the Jacobian
%
@@ -11,12 +11,8 @@ function []=display_problematic_vars_Jacobian(problemrow,problemcol,M_,x,type,ca
% type [string] 'static' or 'dynamic' depending on the type of
% Jacobian
% caller_string [string] contains name of calling function for printing
-%
-% OUTPUTS
-% none.
-%
-% Copyright © 2014-2023 Dynare Team
+% Copyright © 2014-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -40,20 +36,24 @@ end
initial_aux_eq_nbr=M_.ramsey_orig_endo_nbr;
if strcmp(type,'dynamic')
for ii=1:length(problemrow)
- if problemcol(ii)>max(M_.lead_lag_incidence)
- var_row=2;
- var_index=problemcol(ii)-max(max(M_.lead_lag_incidence));
- else
- [var_row,var_index]=find(M_.lead_lag_incidence==problemcol(ii));
- end
- if var_row==2
- type_string='';
- elseif var_row==1
- type_string='lag of';
- elseif var_row==3
- type_string='lead of';
+ if problemcol(ii) <= 3*M_.endo_nbr % Endogenous
+ endo = true;
+ var_index = mod(problemcol(ii)-1, M_.endo_nbr)+1;
+ if problemcol(ii) <= M_.endo_nbr
+ type_string='lag of';
+ elseif problemcol(ii) <= 2*M_.endo_nbr
+ type_string='';
+ else
+ type_string='lead of';
+ end
+ elseif problemcol(ii) <= 3*M_.endo_nbr+M_.exo_nbr % Exogenous
+ endo = false;
+ var_index = problemcol(ii) - 3*M_.endo_nbr;
+ else % Exogenous deterministic, ignore
+ continue
end
- if problemcol(ii)<=max(max(M_.lead_lag_incidence)) && var_index<=M_.orig_endo_nbr
+
+ if endo && var_index <= M_.orig_endo_nbr
if problemrow(ii)<=initial_aux_eq_nbr
eq_nbr = problemrow(ii);
fprintf('Derivative of Auxiliary Equation %d with respect to %s Variable %s (initial value of %s: %g) \n', ...
@@ -63,16 +63,16 @@ if strcmp(type,'dynamic')
fprintf('Derivative of Equation %d with respect to %s Variable %s (initial value of %s: %g) \n', ...
eq_nbr, type_string, M_.endo_names{var_index}, M_.endo_names{var_index}, x(var_index));
end
- elseif problemcol(ii)<=max(max(M_.lead_lag_incidence)) && var_index>M_.orig_endo_nbr % auxiliary vars
- if M_.aux_vars(1,problemcol(ii)-M_.orig_endo_nbr).type==6 %Ramsey Lagrange Multiplier
+ elseif endo && var_index > M_.orig_endo_nbr % auxiliary vars
+ if M_.aux_vars(1, var_index - M_.orig_endo_nbr).type==6 %Ramsey Lagrange Multiplier
if problemrow(ii)<=initial_aux_eq_nbr
eq_nbr = problemrow(ii);
- fprintf('Derivative of Auxiliary Equation %d with respect to %s of Langrange multiplier of equation %s (initial value: %g) \n', ...
- eq_nbr, type_string, M_.aux_vars(1,problemcol(ii)-M_.orig_endo_nbr).eq_nbr, x(problemcol(ii)));
+ fprintf('Derivative of Auxiliary Equation %d with respect to %s of Lagrange multiplier of equation %s (initial value: %g) \n', ...
+ eq_nbr, type_string, M_.aux_vars(1, var_index-M_.orig_endo_nbr).eq_nbr, x(var_index));
else
eq_nbr = problemrow(ii)-initial_aux_eq_nbr;
- fprintf('Derivative of Equation %d with respect to %s of Langrange multiplier of equation %s (initial value: %g) \n', ...
- eq_nbr, type_string, M_.aux_vars(1,problemcol(ii)-M_.orig_endo_nbr).eq_nbr, x(problemcol(ii)));
+ fprintf('Derivative of Equation %d with respect to %s of Lagrange multiplier of equation %s (initial value: %g) \n', ...
+ eq_nbr, type_string, M_.aux_vars(1, var_index-M_.orig_endo_nbr).eq_nbr, x(var_index));
end
else
if problemrow(ii)<=initial_aux_eq_nbr
@@ -87,7 +87,7 @@ if strcmp(type,'dynamic')
eq_nbr, type_string, M_.endo_names{orig_var_index}, M_.endo_names{orig_var_index}, x(orig_var_index));
end
end
- elseif problemcol(ii)>max(max(M_.lead_lag_incidence)) && var_index<=M_.exo_nbr
+ else % Exogenous
if problemrow(ii)<=initial_aux_eq_nbr
eq_nbr = problemrow(ii);
fprintf('Derivative of Auxiliary Equation %d with respect to %s shock %s \n', ...
@@ -97,8 +97,6 @@ if strcmp(type,'dynamic')
fprintf('Derivative of Equation %d with respect to %s shock %s \n', ...
eq_nbr, type_string, M_.exo_names{var_index});
end
- else
- error('display_problematic_vars_Jacobian:: The error should not happen. Please contact the developers')
end
end
fprintf('\n%s The problem most often occurs, because a variable with\n', caller_string)
@@ -149,4 +147,4 @@ elseif strcmp(type, 'static')
fprintf('%s If you are using model-local variables (# operator), check their values as well.\n', caller_string)
else
error('Unknown Type')
-end
\ No newline at end of file
+end
diff --git a/matlab/model_diagnostics.m b/matlab/model_diagnostics.m
index a01947939e..66858d7600 100644
--- a/matlab/model_diagnostics.m
+++ b/matlab/model_diagnostics.m
@@ -16,7 +16,7 @@ function model_diagnostics(M_,options_,oo_)
% none.
%
-% Copyright © 1996-2023 Dynare Team
+% Copyright © 1996-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -271,51 +271,53 @@ if singularity_problem
end
%%check dynamic Jacobian
-klen = M_.maximum_lag + M_.maximum_lead + 1;
-exo_simul = [repmat(oo_.exo_steady_state',klen,1) repmat(oo_.exo_det_steady_state',klen,1)];
-iyv = M_.lead_lag_incidence';
-iyv = iyv(:);
-iyr0 = find(iyv) ;
-it_ = M_.maximum_lag + 1;
-z = repmat(dr.ys,1,klen);
+dyn_endo_ss = repmat(dr.ys, 3, 1);
if options_.order == 1
if (options_.bytecode)
[~, loc_dr] = bytecode('dynamic','evaluate', M_, options_, z, exo_simul, ...
M_.params, dr.ys, 1);
- jacobia_ = [loc_dr.g1 loc_dr.g1_x loc_dr.g1_xd];
+ % TODO: simplify the following once bytecode MEX has been updated to sparse format
+ g1 = zeros(M_.endo_nbr, 3*M_.endo_nbr+M_.exo_nbr+M_.exo_det_nbr);
+ if M_.maximum_endo_lag > 0
+ g1(:, find(M_.lead_lag_incidence(M_.maximum_endo_lag, :))) = loc_dr.g1(:, 1:M_.nspred);
+ end
+ [~,icurr] = find(M_.lead_lag_incidence(M_.maximum_endo_lag+1, :));
+ g1(:, M_.endo_nbr + icurr) = loc_dr.g1(:, M_.nspred+(1:length(icurr)));
+ if M_.maximum_endo_lead > 0
+ g1(:, 2*M_.endo_nbr + find(M_.lead_lag_incidence(M_.maximum_endo_lag+2, :))) = loc_dr.g1(:, M_.nspred+M_.endo_nbr+(1:M_.nsfwrd));
+ end
+ g1(:, 3*M_.endo_nbr+(1:M_.exo_nbr)) = loc_dr.g1_x;
+ g1(:, 3*M_.endo_nbr+M_.exo_nbr+(1:M_.exo_det_nbr)) = loc_dr.g1_xd;
+ g1 = sparse(g1);
else
- [~,jacobia_] = feval([M_.fname '.dynamic'],z(iyr0),exo_simul, ...
- M_.params, dr.ys, it_);
+ g1 = feval([M_.fname '.sparse.dynamic_g1'], dyn_endo_ss, exo, M_.params, dr.ys, ...
+ M_.dynamic_g1_sparse_rowval, M_.dynamic_g1_sparse_colval, ...
+ M_.dynamic_g1_sparse_colptr);
end
elseif options_.order >= 2
- if (options_.bytecode)
- [~, loc_dr] = bytecode('dynamic','evaluate', M_, options_, z, exo_simul, ...
- M_.params, dr.ys, 1);
- jacobia_ = [loc_dr.g1 loc_dr.g1_x];
- else
- [~,jacobia_,hessian1] = feval([M_.fname '.dynamic'],z(iyr0),...
- exo_simul, ...
- M_.params, dr.ys, it_);
- end
+ [g1, T_order, T] = feval([M_.fname '.sparse.dynamic_g1'], dyn_endo_ss, exo, M_.params, ...
+ dr.ys, M_.dynamic_g1_sparse_rowval, M_.dynamic_g1_sparse_colval, ...
+ M_.dynamic_g1_sparse_colptr);
+ g2_v = feval([M_.fname '.sparse.dynamic_g2'], dyn_endo_ss, exo, M_.params, dr.ys, T_order, T);
end
-if any(any(isinf(jacobia_) | isnan(jacobia_)))
+if any(any(isinf(g1) | isnan(g1)))
problem_dummy=1;
- [infrow,infcol]=find(isinf(jacobia_) | isnan(jacobia_));
+ [infrow,infcol]=find(isinf(g1) | isnan(g1));
fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the dynamic model contains Inf or NaN. The problem arises from: \n\n')
display_problematic_vars_Jacobian(infrow,infcol,M_,dr.ys,'dynamic','MODEL_DIAGNOSTICS: ')
end
-if any(any(~isreal(jacobia_)))
- [imagrow,imagcol]=find(abs(imag(jacobia_))>1e-15);
+if any(any(~isreal(g1)))
+ [imagrow,imagcol]=find(abs(imag(g1))>1e-15);
if ~isempty(imagrow)
problem_dummy=1;
fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the dynamic model contains imaginary parts. The problem arises from: \n\n')
display_problematic_vars_Jacobian(imagrow,imagcol,M_,dr.ys,'dynamic','MODEL_DIAGNOSTICS: ')
end
end
-if exist('hessian1','var')
- if any(any(isinf(hessian1) | isnan(hessian1)))
+if exist('g2_v','var')
+ if any(any(isinf(g2_v) | isnan(g2_v)))
problem_dummy=1;
fprintf('\nMODEL_DIAGNOSTICS: The Hessian of the dynamic model contains Inf or NaN.\n')
end
diff --git a/matlab/stochastic_solver/dyn_first_order_solver.m b/matlab/stochastic_solver/dyn_first_order_solver.m
index 3ec2d2b6b3..7d39307bd9 100644
--- a/matlab/stochastic_solver/dyn_first_order_solver.m
+++ b/matlab/stochastic_solver/dyn_first_order_solver.m
@@ -1,9 +1,9 @@
-function [dr, info] = dyn_first_order_solver(jacobia, M_, dr, options_, task)
-% [dr, info] = dyn_first_order_solver(jacobia, M_, dr, options_, task)
+function [dr, info] = dyn_first_order_solver(g1, M_, dr, options_, task)
+% [dr, info] = dyn_first_order_solver(g1, M_, dr, options_, task)
% Computes the first order reduced form of a DSGE model.
%
% INPUTS
-% - jacobia [double] matrix, the jacobian of the dynamic model.
+% - g1 [double] sparse Jacobian of the dynamic model
% - M_ [struct] Matlab's structre describing the model
% - dr [struct] Matlab's structure describing the reduced form model.
% - options_ [struct] Matlab's structure containing the current state of the options
@@ -63,6 +63,7 @@ if isempty(reorder_jacobian_columns)
npred = M_.npred;
nstatic = M_.nstatic;
ndynamic = M_.ndynamic;
+ nspred = M_.nspred;
nsfwrd = M_.nsfwrd;
k1 = 1:(npred+nboth);
@@ -118,8 +119,10 @@ if isempty(reorder_jacobian_columns)
[~,cols_b] = find(lead_lag_incidence(maximum_lag+1, order_var));
- reorder_jacobian_columns = [nonzeros(lead_lag_incidence(:,order_var)'); ...
- nz+(1:exo_nbr)'];
+ reorder_jacobian_columns = [order_var(nstatic+(1:nspred));
+ M_.endo_nbr+order_var(cols_b); ...
+ 2*M_.endo_nbr+order_var(nstatic+npred+(1:nsfwrd));
+ 3*M_.endo_nbr+(1:exo_nbr)'];
end
dr.ghx = [];
@@ -127,7 +130,7 @@ dr.ghu = [];
dr.state_var = state_var;
-jacobia = jacobia(:,reorder_jacobian_columns);
+jacobia = full(g1(:,reorder_jacobian_columns));
if nstatic > 0
[Q, ~] = qr(jacobia(:,index_s));
diff --git a/matlab/stochastic_solver/dyn_second_order_solver.m b/matlab/stochastic_solver/dyn_second_order_solver.m
index 9a44564286..18c62f6f50 100644
--- a/matlab/stochastic_solver/dyn_second_order_solver.m
+++ b/matlab/stochastic_solver/dyn_second_order_solver.m
@@ -1,7 +1,7 @@
-function dr = dyn_second_order_solver(jacobia,hessian_mat,dr,M_,threads_BC)
+function dr = dyn_second_order_solver(g1, g2, 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})
+%! @deftypefn {Function File} {@var{dr} =} dyn_second_order_solver (@var{g1}, @var{g2}, @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
@@ -15,10 +15,10 @@ function dr = dyn_second_order_solver(jacobia,hessian_mat,dr,M_,threads_BC)
%! @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 g1
+%! Sparse matrix containing the Jacobian of the dynamic model
+%! @item g2
+%! Sparse matrix containing the Hessian of the dynamic model
%! @item dr
%! Matlab's structure describing the reduced form solution of the model.
%! @item M_
@@ -36,7 +36,7 @@ function dr = dyn_second_order_solver(jacobia,hessian_mat,dr,M_,threads_BC)
%! @end deftypefn
%@eod:
-% Copyright © 2001-2020 Dynare Team
+% Copyright © 2001-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -58,22 +58,25 @@ 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);
+[~,cols_b] = find(M_.lead_lag_incidence(2, dr.order_var));
+icurr = M_.endo_nbr + dr.order_var(cols_b);
+ilead = 2*M_.endo_nbr + dr.order_var(M_.nstatic+M_.npred+(1:M_.nsfwrd));
+
+kk1 = [dr.order_var(M_.nstatic+(1:M_.nspred)); icurr; ilead;
+ 3*M_.endo_nbr+(1:M_.exo_nbr+M_.exo_det_nbr)'];
+nk = size(g1, 2);
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
+ic = [ M_.nstatic+(1:M_.nspred) ]';
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,:);
+A(:,kcurr~=0) = g1(:, icurr);
+A(:,ic) = A(:,ic) + g1(:, ilead)*dr.ghx(klead~=0,:);
B = zeros(M_.endo_nbr,M_.endo_nbr);
-B(:,M_.nstatic+M_.npred+1:end) = jacobia(:,nonzeros(klead));
+B(:,M_.nstatic+M_.npred+1:end) = g1(:, ilead);
C = dr.ghx(ic,:);
zx = [eye(length(ic));
dr.ghx(kcurr~=0,:);
@@ -85,14 +88,14 @@ zu = [zeros(length(ic),M_.exo_nbr);
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 = sparse_hessian_times_B_kronecker_C(g2(:,kk2(kk1,kk1)), zx, threads_BC); % g2: 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
+rhs = sparse_hessian_times_B_kronecker_C(g2(:,kk2(kk1,kk1)), zx, zu, threads_BC); % g2: reordering to DR order
abcOut = A_times_B_kronecker_C(dr.ghxx, dr.ghx(ic,:), dr.ghu(ic,:));
rhs = -rhs-B*abcOut;
%lhs
@@ -100,7 +103,7 @@ 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
+rhs = sparse_hessian_times_B_kronecker_C(g2(:,kk2(kk1,kk1)), zu, threads_BC); % g2: reordering to DR order
B1 = A_times_B_kronecker_C(B*dr.ghxx,dr.ghu(ic,:));
rhs = -rhs-B1;
%lhs
@@ -111,13 +114,13 @@ dr.ghuu = A\rhs;
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));
+LHS(:,kcurr~=0) = g1(:, icurr);
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;
+B1 = sparse_hessian_times_B_kronecker_C(g2(:,kk2(ilead,ilead)), dr.ghu(klead~=0,:), threads_BC); % g2: focus only on forward variables and reorder to DR order
+RHS = RHS + g1(:, ilead)*dr.ghuu(klead~=0,:)+B1;
% LHS
-LHS = LHS + jacobia(:,nonzeros(klead))*(E(klead~=0,:)+[O1(klead~=0,:) dr.ghx(klead~=0,:) O2(klead~=0,:)]);
+LHS = LHS + g1(:, ilead)*(E(klead~=0,:)+[O1(klead~=0,:) dr.ghx(klead~=0,:) O2(klead~=0,:)]);
RHS = RHS*M_.Sigma_e(:);
dr.fuu = RHS;
RHS = -RHS;
diff --git a/matlab/stochastic_solver/stochastic_solvers.m b/matlab/stochastic_solver/stochastic_solvers.m
index 678810d56e..d382d34e45 100644
--- a/matlab/stochastic_solver/stochastic_solvers.m
+++ b/matlab/stochastic_solver/stochastic_solvers.m
@@ -96,96 +96,100 @@ if options_.k_order_solver
return
end
-klen = M_.maximum_lag + M_.maximum_lead + 1;
-exo_simul = [repmat(exo_steady_state',klen,1) repmat(exo_det_steady_state',klen,1)];
-iyv = M_.lead_lag_incidence';
-iyv = iyv(:);
-iyr0 = find(iyv) ;
+dyn_endo_ss = repmat(dr.ys, 3, 1);
+exo_ss = [exo_steady_state; exo_det_steady_state];
-it_ = M_.maximum_lag + 1;
-z = repmat(dr.ys,1,klen);
if local_order == 1
if (options_.bytecode)
+ klen = M_.maximum_lag + M_.maximum_lead + 1;
+ exo_simul = repmat(exo_ss', klen, 1);
+ z = repmat(dr.ys, 1, klen);
+
[~, loc_dr] = bytecode('dynamic','evaluate', M_, options_, z, exo_simul, ...
M_.params, dr.ys, 1);
- jacobia_ = [loc_dr.g1 loc_dr.g1_x loc_dr.g1_xd];
+ % TODO: simplify the following once bytecode MEX has been updated to sparse format
+ g1 = zeros(M_.endo_nbr, 3*M_.endo_nbr+M_.exo_nbr+M_.exo_det_nbr);
+ if M_.maximum_endo_lag > 0
+ g1(:, find(M_.lead_lag_incidence(M_.maximum_endo_lag, :))) = loc_dr.g1(:, 1:M_.nspred);
+ end
+ [~,icurr] = find(M_.lead_lag_incidence(M_.maximum_endo_lag+1, :));
+ g1(:, M_.endo_nbr + icurr) = loc_dr.g1(:, M_.nspred+(1:length(icurr)));
+ if M_.maximum_endo_lead > 0
+ g1(:, 2*M_.endo_nbr + find(M_.lead_lag_incidence(M_.maximum_endo_lag+2, :))) = loc_dr.g1(:, M_.nspred+M_.endo_nbr+(1:M_.nsfwrd));
+ end
+ g1(:, 3*M_.endo_nbr+(1:M_.exo_nbr)) = loc_dr.g1_x;
+ g1(:, 3*M_.endo_nbr+M_.exo_nbr+(1:M_.exo_det_nbr)) = loc_dr.g1_xd;
+ g1 = sparse(g1);
else
- [~,jacobia_] = feval([M_.fname '.dynamic'],z(iyr0),exo_simul, ...
- M_.params, dr.ys, it_);
+ g1 = feval([M_.fname '.sparse.dynamic_g1'], dyn_endo_ss, exo_ss, M_.params, dr.ys, ...
+ M_.dynamic_g1_sparse_rowval, M_.dynamic_g1_sparse_colval, ...
+ M_.dynamic_g1_sparse_colptr);
end
elseif local_order == 2
if (options_.bytecode)
warning('Option "bytecode" is ignored when computing perturbation solution at order = 2')
end
- [~,jacobia_,hessian1] = feval([M_.fname '.dynamic'],z(iyr0),...
- exo_simul, ...
- M_.params, dr.ys, it_);
- [infrow, ~] = find(isinf(hessian1));
- if options_.debug
- if ~isempty(infrow)
+ [g1, T_order, T] = feval([M_.fname '.sparse.dynamic_g1'], dyn_endo_ss, exo_ss, M_.params, ...
+ dr.ys, M_.dynamic_g1_sparse_rowval, M_.dynamic_g1_sparse_colval, ...
+ M_.dynamic_g1_sparse_colptr);
+ g2_v = feval([M_.fname '.sparse.dynamic_g2'], dyn_endo_ss, exo_ss, M_.params, dr.ys, T_order, T);
+
+ g2 = build_two_dim_hessian(M_.dynamic_g2_sparse_indices, g2_v, size(g1, 1), size(g1, 2));
+
+ if any(any(isinf(g2)))
+ if options_.debug
fprintf('\nSTOCHASTIC_SOLVER: The Hessian of the dynamic model contains Inf.\n')
fprintf('STOCHASTIC_SOLVER: Try running model_diagnostics to find the source of the problem.\n')
- save([M_.dname filesep 'Output' filesep M_.fname '_debug.mat'],'hessian1')
+ save([M_.dname filesep 'Output' filesep M_.fname '_debug.mat'], 'g2')
end
- end
- if ~isempty(infrow)
info(1)=11;
return
end
- [nanrow, ~] = find(isnan(hessian1));
- if options_.debug
- if ~isempty(nanrow)
+
+ if any(any(isnan(g2)))
+ if options_.debug
fprintf('\nSTOCHASTIC_SOLVER: The Hessian of the dynamic model contains NaN.\n')
fprintf('STOCHASTIC_SOLVER: Try running model_diagnostics to find the source of the problem.\n')
- save([M_.dname filesep 'Output' filesep M_.fname '_debug.mat'],'hessian1')
+ save([M_.dname filesep 'Output' filesep M_.fname '_debug.mat'], 'g2')
end
- end
- if ~isempty(nanrow)
info(1)=12;
return
end
end
-[infrow,infcol]=find(isinf(jacobia_));
-
-if options_.debug
- if ~isempty(infrow)
+[infrow, infcol] = find(isinf(g1));
+if ~isempty(infrow)
+ if options_.debug
fprintf('\nSTOCHASTIC_SOLVER: The Jacobian of the dynamic model contains Inf. The problem is associated with:\n\n')
display_problematic_vars_Jacobian(infrow,infcol,M_,dr.ys,'dynamic','STOCHASTIC_SOLVER: ')
- save([M_.dname filesep 'Output' filesep M_.fname '_debug.mat'],'jacobia_')
+ save([M_.dname filesep 'Output' filesep M_.fname '_debug.mat'], 'g1')
end
-end
-
-if ~isempty(infrow)
info(1)=10;
return
end
-if ~isreal(jacobia_)
- if max(max(abs(imag(jacobia_)))) < 1e-15
- jacobia_ = real(jacobia_);
+if ~isreal(g1)
+ if max(max(abs(imag(g1)))) < 1e-15
+ g1 = real(g1);
else
if options_.debug
- [imagrow,imagcol]=find(abs(imag(jacobia_))>1e-15);
+ [imagrow, imagcol]=find(abs(imag(g1)) > 1e-15);
fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the dynamic model contains imaginary parts. The problem arises from: \n\n')
display_problematic_vars_Jacobian(imagrow,imagcol,M_,dr.ys,'dynamic','STOCHASTIC_SOLVER: ')
end
info(1) = 6;
- info(2) = sum(sum(imag(jacobia_).^2));
+ info(2) = sum(sum(imag(g1).^2));
return
end
end
-[nanrow,nancol]=find(isnan(jacobia_));
-if options_.debug
- if ~isempty(nanrow)
+[nanrow, nancol] = find(isnan(g1));
+if ~isempty(nanrow)
+ if options_.debug
fprintf('\nSTOCHASTIC_SOLVER: The Jacobian of the dynamic model contains NaN. The problem is associated with:\n\n')
display_problematic_vars_Jacobian(nanrow,nancol,M_,dr.ys,'dynamic','STOCHASTIC_SOLVER: ')
- save([M_.dname filesep 'Output' filesep M_.fname '_debug.mat'],'jacobia_')
+ save([M_.dname filesep 'Output' filesep M_.fname '_debug.mat'], 'g1')
end
-end
-
-if ~isempty(nanrow)
info(1) = 8;
NaN_params=find(isnan(M_.params));
info(2:length(NaN_params)+1) = NaN_params;
@@ -193,20 +197,14 @@ if ~isempty(nanrow)
end
nstatic = M_.nstatic;
+npred = M_.npred;
nfwrd = M_.nfwrd;
nspred = M_.nspred;
nboth = M_.nboth;
nsfwrd = M_.nsfwrd;
order_var = dr.order_var;
nd = M_.nspred+M_.nsfwrd;
-nz = nnz(M_.lead_lag_incidence);
-
-sdyn = M_.endo_nbr - nstatic;
-
-[~,cols_b,cols_j] = find(M_.lead_lag_incidence(M_.maximum_endo_lag+1, ...
- order_var));
-b = zeros(M_.endo_nbr,M_.endo_nbr);
-b(:,cols_b) = jacobia_(:,cols_j);
+[~,icurr_dr] = find(M_.lead_lag_incidence(M_.maximum_endo_lag+1, order_var));
if M_.maximum_endo_lead == 0
% backward models: simplified code exist only at order == 1
@@ -216,9 +214,10 @@ if M_.maximum_endo_lead == 0
else
dr.state_var = [];
end
- dr.ghx = -b\jacobia_(:,1:nspred);
+ b = full(g1(:, M_.endo_nbr+order_var(icurr_dr)));
+ dr.ghx = -b \ full(g1(:, order_var(nstatic+(1:nspred))));
if M_.exo_nbr
- dr.ghu = -b\jacobia_(:,nz+1:end);
+ dr.ghu = -b \ full(g1(:, 3*M_.endo_nbr+1:end));
end
dr.eigval = eig(kalman_transition_matrix(dr,nstatic+(1:nspred),1:nspred));
dr.full_rank = 1;
@@ -240,10 +239,9 @@ if M_.maximum_endo_lead == 0
else
% If required, use AIM solver if not check only
if options_.aim_solver && (task == 0)
- [dr,info] = AIM_first_order_solver(jacobia_,M_,dr,options_.qz_criterium);
-
+ [dr, info] = AIM_first_order_solver(g1, M_, dr, options_.qz_criterium);
else % use original Dynare solver
- [dr,info] = dyn_first_order_solver(jacobia_,M_,dr,options_,task);
+ [dr, info] = dyn_first_order_solver(g1, M_, dr, options_, task);
if info(1) || task
return
end
@@ -251,17 +249,8 @@ else
if local_order > 1
% Second order
- dr = dyn_second_order_solver(jacobia_,hessian1,dr,M_,...
+ dr = dyn_second_order_solver(g1, g2, dr, M_, ...
options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
-
- % reordering second order derivatives, used for deterministic
- % variables below
- k1 = nonzeros(M_.lead_lag_incidence(:,order_var)');
- kk = [k1; length(k1)+(1:M_.exo_nbr+M_.exo_det_nbr)'];
- nk = size(kk,1);
- kk1 = reshape(1:nk^2,nk,nk);
- kk1 = kk1(kk,kk);
- hessian1 = hessian1(:,kk1(:));
end
end
@@ -269,9 +258,9 @@ end
%exogenous deterministic variables
if M_.exo_det_nbr > 0
gx = dr.gx;
- f1 = sparse(jacobia_(:,nonzeros(M_.lead_lag_incidence(M_.maximum_endo_lag+2:end,order_var))));
- f0 = sparse(jacobia_(:,nonzeros(M_.lead_lag_incidence(M_.maximum_endo_lag+1,order_var))));
- fudet = sparse(jacobia_(:,nz+M_.exo_nbr+1:end));
+ f1 = g1(:,2*M_.endo_nbr + order_var(nstatic+npred+(1:nsfwrd)));
+ f0 = g1(:,M_.endo_nbr + order_var(icurr_dr));
+ fudet = g1(:,3*M_.endo_nbr+M_.exo_nbr+1:end);
M1 = inv(f0+[zeros(M_.endo_nbr,nstatic) f1*gx zeros(M_.endo_nbr,nsfwrd-nboth)]);
M2 = M1*f1;
dr.ghud = cell(M_.exo_det_length,1);
@@ -281,6 +270,15 @@ if M_.exo_det_nbr > 0
end
if local_order > 1
+ % reordering second order derivatives
+ kk1 = [order_var(nstatic+(1:nspred));
+ M_.endo_nbr + order_var(icurr_dr);
+ 2*M_.endo_nbr + order_var(nstatic+npred+(1:nsfwrd));
+ 3*M_.endo_nbr+(1:M_.exo_nbr+M_.exo_det_nbr)'];
+ nk = size(g1, 2);
+ kk2 = reshape(1:nk^2, nk, nk);
+ g2_reordered = g2(:,kk2(kk1,kk1));
+
lead_lag_incidence = M_.lead_lag_incidence;
k0 = find(lead_lag_incidence(M_.maximum_endo_lag+1,order_var)');
k1 = find(lead_lag_incidence(M_.maximum_endo_lag+2,order_var)');
@@ -291,7 +289,7 @@ if M_.exo_det_nbr > 0
zu = [zeros(nspred,M_.exo_nbr); dr.ghu(k0,:); gx*hu; zeros(M_.exo_nbr+M_.exo_det_nbr, ...
M_.exo_nbr)];
zud=[zeros(nspred,M_.exo_det_nbr);dr.ghud{1};gx(:,1:nspred)*hud;zeros(M_.exo_nbr,M_.exo_det_nbr);eye(M_.exo_det_nbr)];
- R1 = hessian1*kron(zx,zud);
+ R1 = g2_reordered*kron(zx,zud);
dr.ghxud = cell(M_.exo_det_length,1);
kf = M_.endo_nbr-nfwrd-nboth+1:M_.endo_nbr;
kp = nstatic+[1:nspred];
@@ -300,37 +298,37 @@ if M_.exo_det_nbr > 0
for i = 2:M_.exo_det_length
hudi = dr.ghud{i}(kp,:);
zudi=[zeros(nspred,M_.exo_det_nbr);dr.ghud{i};gx(:,1:nspred)*hudi;zeros(M_.exo_nbr+M_.exo_det_nbr,M_.exo_det_nbr)];
- R2 = hessian1*kron(zx,zudi);
+ R2 = g2_reordered*kron(zx,zudi);
dr.ghxud{i} = -M2*(dr.ghxud{i-1}(kf,:)*kron(dr.Gy,Eud)+dr.ghxx(kf,:)*kron(dr.ghx(kp,:),dr.ghud{i}(kp,:)))-M1*R2;
end
- R1 = hessian1*kron(zu,zud);
+ R1 = g2_reordered*kron(zu,zud);
dr.ghudud = cell(M_.exo_det_length,1);
dr.ghuud{1} = -M1*(R1+f1*dr.ghxx(kf,:)*kron(dr.ghu(kp,:),dr.ghud{1}(kp,:)));
Eud = eye(M_.exo_det_nbr);
for i = 2:M_.exo_det_length
hudi = dr.ghud{i}(kp,:);
zudi=[zeros(nspred,M_.exo_det_nbr);dr.ghud{i};gx(:,1:nspred)*hudi;zeros(M_.exo_nbr+M_.exo_det_nbr,M_.exo_det_nbr)];
- R2 = hessian1*kron(zu,zudi);
+ R2 = g2_reordered*kron(zu,zudi);
dr.ghuud{i} = -M2*dr.ghxud{i-1}(kf,:)*kron(hu,Eud)-M1*R2;
end
- R1 = hessian1*kron(zud,zud);
+ R1 = g2_reordered*kron(zud,zud);
dr.ghudud = cell(M_.exo_det_length,M_.exo_det_length);
dr.ghudud{1,1} = -M1*R1-M2*dr.ghxx(kf,:)*kron(hud,hud);
for i = 2:M_.exo_det_length
hudi = dr.ghud{i}(nstatic+1:nstatic+nspred,:);
zudi=[zeros(nspred,M_.exo_det_nbr);dr.ghud{i};gx(:,1:nspred)*hudi+dr.ghud{i-1}(kf,:);zeros(M_.exo_nbr+M_.exo_det_nbr,M_.exo_det_nbr)];
- R2 = hessian1*kron(zudi,zudi);
+ R2 = g2_reordered*kron(zudi,zudi);
dr.ghudud{i,i} = -M2*(dr.ghudud{i-1,i-1}(kf,:)+...
2*dr.ghxud{i-1}(kf,:)*kron(hudi,Eud) ...
+dr.ghxx(kf,:)*kron(hudi,hudi))-M1*R2;
- R2 = hessian1*kron(zud,zudi);
+ R2 = g2_reordered*kron(zud,zudi);
dr.ghudud{1,i} = -M2*(dr.ghxud{i-1}(kf,:)*kron(hud,Eud)+...
dr.ghxx(kf,:)*kron(hud,hudi))...
-M1*R2;
for j=2:i-1
hudj = dr.ghud{j}(kp,:);
zudj=[zeros(nspred,M_.exo_det_nbr);dr.ghud{j};gx(:,1:nspred)*hudj;zeros(M_.exo_nbr+M_.exo_det_nbr,M_.exo_det_nbr)];
- R2 = hessian1*kron(zudj,zudi);
+ R2 = g2_reordered*kron(zudj,zudi);
dr.ghudud{j,i} = -M2*(dr.ghudud{j-1,i-1}(kf,:)+dr.ghxud{j-1}(kf,:)* ...
kron(hudi,Eud)+dr.ghxud{i-1}(kf,:)* ...
kron(hudj,Eud)+dr.ghxx(kf,:)*kron(hudj,hudi))-M1*R2;
--
GitLab