From e8c4c191cb31cb4d9666ec63ba37d821de456ccb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Tue, 9 Apr 2024 16:46:10 +0200
Subject: [PATCH] Consolidate AIM functions
---
license.txt | 2 +-
matlab/AIM/dynAIMsolver1.m | 111 ------------------
matlab/get_error_message.m | 7 +-
.../AIM_first_order_solver.m | 101 +++++++++++++++-
.../stochastic_solver/convertAimCodeToInfo.m | 64 ----------
5 files changed, 101 insertions(+), 184 deletions(-)
delete mode 100644 matlab/AIM/dynAIMsolver1.m
delete mode 100644 matlab/stochastic_solver/convertAimCodeToInfo.m
diff --git a/license.txt b/license.txt
index c16b157650..45842121d5 100644
--- a/license.txt
+++ b/license.txt
@@ -64,7 +64,7 @@ License: public-domain-occbin
dynamic models with occasionally binding constraints easily"
Journal of Monetary Economics 70, 22-38
-Files: matlab/AIM/SP*
+Files: matlab/AIM/*
Copyright: none
License: public-domain-aim
This code is in the public domain and may be used freely.
diff --git a/matlab/AIM/dynAIMsolver1.m b/matlab/AIM/dynAIMsolver1.m
deleted file mode 100644
index 7c02fbf531..0000000000
--- a/matlab/AIM/dynAIMsolver1.m
+++ /dev/null
@@ -1,111 +0,0 @@
-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, . . . ,?
-% and its input as single array of matrices: [H-$... Hi ... H+&]
-% and its solution as xt=SUM( Bi*xt+i) + @*£*zt for i=-$...-1
-% with the output in form bb=[B-$... Bi ... B-1] and @=inv(Ho+H1*B-1)
-% Dynare jacobian = [fy'-$... fy'i ... fy'+& fu']
-% where [fy'-$... fy'i ... fy'+&]=[H-$... Hi ... H+&] and fu'= £
-%
-% INPUTS
-% g1 [matrix] sparse Jacobian of the dynamic model
-% dr [matlab structure] Decision rules for stochastic simulations.
-% M_ [matlab structure] Definition of the model.
-%
-% OUTPUTS
-% dr [matlab structure] Decision rules for stochastic simulations.
-% aimcode [integer] 1: the model defines variables uniquely
-% aimcode is resolved in AIMerr as
-% (c==1) e='Aim: unique solution.';
-% (c==2) e='Aim: roots not correctly computed by real_schur.';
-% (c==3) e='Aim: too many big roots.';
-% (c==35) e='Aim: too many big roots, and q(:,right) is singular.';
-% (c==4) e='Aim: too few big roots.';
-% (c==45) e='Aim: too few big roots, and q(:,right) is singular.';
-% (c==5) e='Aim: q(:,right) is singular.';
-% (c==61) e='Aim: too many exact shiftrights.';
-% (c==62) e='Aim: too many numeric shiftrights.';
-% (c==63) e='Aim: A is NAN or INF.';
-% (c==64) e='Aim: Problem in SPEIG.';
-% else e='Aimerr: return code not properly specified';
-%
-% SPECIAL REQUIREMENTS
-% Dynare use:
-% 1) the lognormal block in DR1 is being invoked for some models and changing
-% values of ghx and ghy. We need to return the AIM output
-% values before that block and run the block with the current returned values
-% of gy (i.e. dr.ghx) and gu (dr.ghu) if it is needed even when the AIM is used
-% (it does not depend on mjdgges output).
-%
-% 2) passing in aa={Q'|1}*jacobia_ can produce ~ one order closer
-% results to the Dynare solutiion then when if plain jacobia_ is passed,
-% i.e. diff < e-14 for aa and diff < *e-13 for jacobia_ if Q' is used.
-%
-% GP July 2008
-
-% Copyright © 2008-2024 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/>.
-
-aimcode=-1;
-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
-
-%disp('gensysToAMA:running ama');
-try % try to run AIM
- [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
- dr.ghx = bb(dr.order_var, dr.order_var(M_.nstatic+(1:M_.nspred)));
-
- 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);
- %? = 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(g1(:, 1:3*M_.endo_nbr), bb, neq, 1, 1);
- scof1 = scof(:, neq+1:end);
- scof1= scof1(:,dr.order_var);
- dr.ghu = -scof1 \ g1(:, 3*M_.endo_nbr+1:end);
- else
- dr.ghu = [];
- end
-else
- err=SPAimerr(aimcode);
- %warning('Error in AIM: aimcode=%d, erro=%s', aimcode, err);;
- disp(['Error in AIM: aimcode=' sprintf('%d : %s',aimcode, err)]);
- if aimcode < 1 || aimcode > 5 % too big exception, use mjdgges
- error('Error in AIM: aimcode=%d ; %s', aimcode, err);
- end
- % if aimcode > 5
- % disp(['Error in AIM: aimcode=' sprintf('%d : %s',aimcode, err)]);
- % aimcode=5;
- % end
-end
diff --git a/matlab/get_error_message.m b/matlab/get_error_message.m
index 1694dde4e2..9680a67573 100644
--- a/matlab/get_error_message.m
+++ b/matlab/get_error_message.m
@@ -10,7 +10,7 @@ function message = get_error_message(info, options_)
% SPECIAL REQUIREMENTS
% none
-% Copyright © 2005-2023 Dynare Team
+% Copyright © 2005-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -129,7 +129,6 @@ switch info(1)
message = 'Calibrated covariance of the structural errors implies correlation larger than +-1.';
case 72
message = 'Calibrated covariance of the measurement errors implies correlation larger than +-1.';
- % Aim Code Conversions by convertAimCodeToInfo.m
case 81
message = ['Ramsey: The solution to the static first order conditions for optimal policy could not be found. Either the model' ...
' doesn''t have a steady state, there are an infinity of steady states, ' ...
@@ -163,7 +162,7 @@ switch info(1)
case 162
message = 'Aim: too many numeric shiftrights';
case 163
- message = 'Aim: A is NAN or INF.';
+ message = 'Aim: A is NaN or Inf.';
case 164
message = 'Aim: Problem in SPEIG.';
case 180
@@ -212,4 +211,4 @@ switch info(1)
message = 'Logarithmic reduction converged to a solution that does not solve the matrix equation.';
otherwise
message = 'This case shouldn''t happen. Contact the authors of Dynare';
-end
\ No newline at end of file
+end
diff --git a/matlab/stochastic_solver/AIM_first_order_solver.m b/matlab/stochastic_solver/AIM_first_order_solver.m
index a8e48fbd7d..067d554a35 100644
--- a/matlab/stochastic_solver/AIM_first_order_solver.m
+++ b/matlab/stochastic_solver/AIM_first_order_solver.m
@@ -1,7 +1,7 @@
-function [dr,info]=AIM_first_order_solver(jacobia,M_,dr,qz_criterium)
+function [dr,info]=AIM_first_order_solver(g1,M_,dr,qz_criterium)
%@info:
-%! @deftypefn {Function File} {[@var{dr},@var{info}] =} AIM_first_order_solver (@var{jacobia},@var{M_},@var{dr},@var{qz_criterium})
+%! @deftypefn {Function File} {[@var{dr},@var{info}] =} AIM_first_order_solver (@var{g1},@var{M_},@var{dr},@var{qz_criterium})
%! @anchor{AIM_first_order_solver}
%! @sp 1
%! Computes the first order reduced form of the DSGE model using AIM.
@@ -9,7 +9,7 @@ function [dr,info]=AIM_first_order_solver(jacobia,M_,dr,qz_criterium)
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
-%! @item jacobia
+%! @item g1
%! Matrix containing the Jacobian of the model
%! @item M_
%! Matlab's structure describing the model (initialized by @code{dynare}).
@@ -50,6 +50,29 @@ function [dr,info]=AIM_first_order_solver(jacobia,M_,dr,qz_criterium)
%! @end table
%! @end deftypefn
%@eod:
+%
+% Maps Dynare jacobian to AIM 1st order model solver designed and developed by Gary Anderson
+% and derives the solution for dr.ghx and dr.ghu from the AIM outputs
+% AIM System is given as a sum:
+% i.e. for i=-$...+& SUM(Hi*xt+i)= £*zt, t = 0, . . . ,?
+% and its input as single array of matrices: [H-$... Hi ... H+&]
+% and its solution as xt=SUM( Bi*xt+i) + @*£*zt for i=-$...-1
+% with the output in form bb=[B-$... Bi ... B-1] and @=inv(Ho+H1*B-1)
+% Dynare jacobian = [fy'-$... fy'i ... fy'+& fu']
+% where [fy'-$... fy'i ... fy'+&]=[H-$... Hi ... H+&] and fu'= £
+%
+% Dynare use:
+% 1) the lognormal block in DR1 is being invoked for some models and changing
+% values of ghx and ghy. We need to return the AIM output
+% values before that block and run the block with the current returned values
+% of dr.ghx and dr.ghu if it is needed even when the AIM is used
+% (it does not depend on mjdgges output).
+%
+% 2) passing in aa={Q'|1}*g1 can produce ~ one order closer
+% results to the Dynare solutiion then when if plain g1 is passed,
+% i.e. diff < e-14 for aa and diff < *e-13 for g1 if Q' is used.
+%
+% Initially written by George Perendia
% Copyright © 2001-2024 Dynare Team
%
@@ -70,7 +93,46 @@ function [dr,info]=AIM_first_order_solver(jacobia,M_,dr,qz_criterium)
info = 0;
-[dr,aimcode]=dynAIMsolver1(jacobia,M_,dr);
+aimcode=-1;
+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
+
+try % try to run AIM
+ [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
+ dr.ghx = bb(dr.order_var, dr.order_var(M_.nstatic+(1:M_.nspred)));
+
+ if M_.exo_nbr % if there are exogenous shocks then derive ghu 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);
+ %? = 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 ghu
+ % Using AIM SPObstruct
+ 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 \ g1(:, 3*M_.endo_nbr+1:end);
+ else
+ dr.ghu = [];
+ end
+else
+ err=SPAimerr(aimcode);
+ disp(['Error in AIM: aimcode=' sprintf('%d : %s',aimcode, err)]);
+ if aimcode < 1 || aimcode > 5 % too big exception, use mjdgges
+ error('Error in AIM: aimcode=%d ; %s', aimcode, err);
+ end
+end
if aimcode ~=1
info(1) = convertAimCodeToInfo(aimCode); %convert to be in the 100 range
@@ -97,3 +159,34 @@ if nba ~= nsfwrd
info(2) = temp'*temp;
return
end
+
+
+function [info] = convertAimCodeToInfo(aimCode)
+% Returns an appropriate code for get_error_message.m (see that function for the meaning)
+
+switch aimCode
+ case 1
+ info = 0;
+ case 2
+ info = 102;
+ case 3
+ info = 103;
+ case 35
+ info = 135;
+ case 4
+ info = 104;
+ case 45
+ info = 145;
+ case 5
+ info = 105;
+ case 61
+ info = 161;
+ case 62
+ info = 162;
+ case 63
+ info = 163;
+ case 64
+ info = 164;
+ otherwise
+ info = 1;
+end
diff --git a/matlab/stochastic_solver/convertAimCodeToInfo.m b/matlab/stochastic_solver/convertAimCodeToInfo.m
deleted file mode 100644
index 243c4c898c..0000000000
--- a/matlab/stochastic_solver/convertAimCodeToInfo.m
+++ /dev/null
@@ -1,64 +0,0 @@
-function [info] = convertAimCodeToInfo(aimCode)
-% function [info] = convertAimCodeToInfo(aimCode)
-% Returns an appropriate code for print_info
-%
-% INPUTS
-% aimCode [integer] code returned by AIM
-% (aimCode==1) e='Aim: unique solution.';
-% (aimCode==2) e='Aim: roots not correctly computed by real_schur.';
-% (aimCode==3) e='Aim: too many big roots.';
-% (aimCode==35) e='Aim: too many big roots, and q(:,right) is singular.';
-% (aimCode==4) e='Aim: too few big roots.';
-% (aimCode==45) e='Aim: too few big roots, and q(:,right) is singular.';
-% (aimCode==5) e='Aim: q(:,right) is singular.';
-% (aimCode==61) e='Aim: too many exact shiftrights.';
-% (aimCode==62) e='Aim: too many numeric shiftrights.';
-% (aimCode==63) e='Aim: A is NAN or INF.';
-% (aimCode==64) e='Aim: Problem in SPEIG.';
-%
-% OUTPUTS
-% info [integer] Code to be used to print error in print_info.m
-
-% Copyright © 2011-2017 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/>.
-
-switch aimCode
- case 1
- info = 0; % no problem encountered
- case 2
- info = 102;
- case 3
- info = 103;
- case 35
- info = 135;
- case 4
- info = 104;
- case 45
- info = 145;
- case 5
- info = 105;
- case 61
- info = 161;
- case 62
- info = 162;
- case 63
- info = 163;
- case 64
- info = 164;
- otherwise
- info = 1;
-end
\ No newline at end of file
--
GitLab