function [dr,aimcode]=dynAIMsolver1(jacobia_,M_,dr)
Maps Dynare jacobia 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
jacobia_ [matrix] 1st order derivative of the 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 20080001 function [dr,aimcode,rts]=dynAIMsolver1(jacobia_,M_,dr) 0002 % function [dr,aimcode]=dynAIMsolver1(jacobia_,M_,dr) 0003 % Maps Dynare jacobia to AIM 1st order model solver designed and developed by Gary ANderson 0004 % and derives the solution for gy=dr.hgx and gu=dr.hgu from the AIM outputs 0005 % AIM System is given as a sum: 0006 % i.e. for i=-$...+& SUM(Hi*xt+i)= �*zt, t = 0, . . . ,? 0007 % and its input as single array of matrices: [H-$... Hi ... H+&] 0008 % and its solution as xt=SUM( Bi*xt+i) + @*�*zt for i=-$...-1 0009 % with the output in form bb=[B-$... Bi ... B-1] and @=inv(Ho+H1*B-1) 0010 % Dynare jacobian = [fy'-$... fy'i ... fy'+& fu'] 0011 % where [fy'-$... fy'i ... fy'+&]=[H-$... Hi ... H+&] and fu'= � 0012 % 0013 % INPUTS 0014 % jacobia_ [matrix] 1st order derivative of the model 0015 % dr [matlab structure] Decision rules for stochastic simulations. 0016 % M_ [matlab structure] Definition of the model. 0017 % 0018 % OUTPUTS 0019 % dr [matlab structure] Decision rules for stochastic simulations. 0020 % aimcode [integer] 1: the model defines variables uniquely 0021 % aimcode is resolved in AIMerr as 0022 % (c==1) e='Aim: unique solution.'; 0023 % (c==2) e='Aim: roots not correctly computed by real_schur.'; 0024 % (c==3) e='Aim: too many big roots.'; 0025 % (c==35) e='Aim: too many big roots, and q(:,right) is singular.'; 0026 % (c==4) e='Aim: too few big roots.'; 0027 % (c==45) e='Aim: too few big roots, and q(:,right) is singular.'; 0028 % (c==5) e='Aim: q(:,right) is singular.'; 0029 % (c==61) e='Aim: too many exact shiftrights.'; 0030 % (c==62) e='Aim: too many numeric shiftrights.'; 0031 % (c==63) e='Aim: A is NAN or INF.'; 0032 % (c==64) e='Aim: Problem in SPEIG.'; 0033 % else e='Aimerr: return code not properly specified'; 0034 % 0035 % SPECIAL REQUIREMENTS 0036 % Dynare use: 0037 % 1) the lognormal block in DR1 is being invoked for some models and changing 0038 % values of ghx and ghy. We need to return the AIM output 0039 % values before that block and run the block with the current returned values 0040 % of gy (i.e. dr.ghx) and gu (dr.ghu) if it is needed even when the AIM is used 0041 % (it does not depend on mjdgges output). 0042 % 0043 % 2) passing in aa={Q'|1}*jacobia_ can produce ~ one order closer 0044 % results to the Dynare solutiion then when if plain jacobia_ is passed, 0045 % i.e. diff < e-14 for aa and diff < *e-13 for jacobia_ if Q' is used. 0046 % 0047 % GP July 2008 0048 0049 % Copyright (C) 2008-2009 Dynare Team 0050 % 0051 % This file is part of Dynare. 0052 % 0053 % Dynare is free software: you can redistribute it and/or modify 0054 % it under the terms of the GNU General Public License as published by 0055 % the Free Software Foundation, either version 3 of the License, or 0056 % (at your option) any later version. 0057 % 0058 % Dynare is distributed in the hope that it will be useful, 0059 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0060 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0061 % GNU General Public License for more details. 0062 % 0063 % You should have received a copy of the GNU General Public License 0064 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0065 0066 aimcode=-1; 0067 neq= size(jacobia_,1); % no of equations 0068 lags=M_.maximum_endo_lag; % no of lags and leads 0069 leads=M_.maximum_endo_lead; 0070 klen=(leads+lags+1); % total lenght 0071 theAIM_H=zeros(neq, neq*klen); % alloc space 0072 lli=M_.lead_lag_incidence'; 0073 % "sparse" the compact jacobia into AIM H aray of matrices 0074 % without exogenous shoc 0075 theAIM_H(:,find(lli(:)))=jacobia_(:,nonzeros(lli(:))); 0076 condn = 1.e-10;%SPAmalg uses this in zero tests 0077 uprbnd = 1 + 1.e-6;%allow unit roots 0078 % forward only models - AIM must have at least 1 lead and 1 lag. 0079 if lags ==0 0080 theAIM_H =[zeros(neq) theAIM_H]; 0081 lags=1; 0082 klen=(leads+lags+1); 0083 end 0084 % backward looking only models 0085 if leads ==0 0086 theAIM_H =[theAIM_H zeros(neq)]; 0087 leads=1; 0088 klen=(leads+lags+1); 0089 end 0090 %disp('gensysToAMA:running ama'); 0091 try % try to run AIM 0092 [bb,rts,ia,nexact,nnumeric,lgroots,aimcode] =... 0093 SPAmalg(theAIM_H,neq, lags,leads,condn,uprbnd); 0094 catch 0095 err = lasterror; 0096 disp(['Dynare AIM Solver error:' sprintf('%s; ID:%s',err.message, err.identifier)]); 0097 rethrow(lasterror); 0098 end 0099 if aimcode==1 %if OK 0100 col_order=[]; 0101 for i =1:lags 0102 col_order=[((i-1)*neq)+dr.order_var' col_order]; 0103 end 0104 bb_ord= bb(dr.order_var,col_order); % derive ordered gy 0105 0106 % variables are present in the state space at the lag at which they 0107 % appear and at all smaller lags. The are ordered from smaller to 0108 % higher lag (reversed order of M_.lead_lag_incidence rows for lagged 0109 % variables) 0110 i_lagged_vars = flipud(cumsum(M_.lead_lag_incidence(1:M_.maximum_lag,dr.order_var),1))'; 0111 0112 dr.ghx= bb_ord(:,find(i_lagged_vars(:))); % get columns reported in 0113 % Dynare solution 0114 if M_.exo_nbr % if there are exogenous shocks then derive gu for the shocks: 0115 % get H0 and H+1=HM 0116 % theH0= theAIM_H (:,M_.maximum_endo_lag*neq+1: (M_.maximum_endo_lag+1)*neq); 0117 %theH0= theAIM_H (:,lags*neq+1: (lags+1)*neq); 0118 % theHP= theAIM_H (:,(M_.maximum_endo_lag+1)*neq+1: (M_.maximum_endo_lag+2)*neq); 0119 %theHP= theAIM_H (:,(lags+1)*neq+1: (lags+2)*neq); 0120 theAIM_Psi= - jacobia_(:, size(nonzeros(lli(:)))+1:end);% 0121 %? = inv(H0 + H1B1) 0122 %phi= (theH0+theHP*sparse(bb(:,(lags-1)*neq+1:end)))\eye( neq); 0123 %AIM_ghu=phi*theAIM_Psi; 0124 %dr.ghu =AIM_ghu(dr.order_var,:); % order gu 0125 % Using AIM SPObstruct 0126 scof = SPObstruct(theAIM_H,bb,neq,lags,leads); 0127 scof1= scof(:,(lags)*neq+1:end); 0128 scof1= scof1(:,dr.order_var); 0129 dr.ghu =scof1\theAIM_Psi; 0130 else 0131 dr.ghu = []; 0132 end 0133 else 0134 err=SPAimerr(aimcode); 0135 %warning('Error in AIM: aimcode=%d, erro=%s', aimcode, err);; 0136 disp(['Error in AIM: aimcode=' sprintf('%d : %s',aimcode, err)]); 0137 if aimcode < 1 || aimcode > 5 % too big exception, use mjdgges 0138 error('Error in AIM: aimcode=%d ; %s', aimcode, err); 0139 end 0140 % if aimcode > 5 0141 % disp(['Error in AIM: aimcode=' sprintf('%d : %s',aimcode, err)]); 0142 % aimcode=5; 0143 % end 0144 end