Skip to content
Snippets Groups Projects
Verified Commit e8c4c191 authored by Sébastien Villemot's avatar Sébastien Villemot
Browse files

Consolidate AIM functions

parent a36287fc
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
license.txt
+ 1
− 1
View file @ e8c4c191
......@@ -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.
......
matlab/AIM/dynAIMsolver1.m deleted 100644 → 0
+ 0
− 111
View file @ a36287fc
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
matlab/get_error_message.m
+ 3
− 4
View file @ e8c4c191
......@@ -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
......
matlab/stochastic_solver/AIM_first_order_solver.m
+ 97
− 4
View file @ e8c4c191
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
matlab/stochastic_solver/convertAimCodeToInfo.m deleted 100644 → 0
+ 0
− 64
View file @ a36287fc
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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment