git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1965 ac1d8469-bf42-47a9-8791-bf33cf982152
george
parent
c9cf49e30b
commit
8225da0e7e
|
|
@ -0,0 +1,75 @@
|
|||
AIM Subsystem cosists of
|
||||
|
||||
A function dynAIMsolver1(jacobia_,M_,dr) which maps Dynare to Gary Anderson AIM package subsystem and derives the solution for gy=dr.hgx and gu=dr.hgu from the AIM outputs (where "1" in the title is for 1st order solver).
|
||||
|
||||
A subset of routines from Gary Anderson AIM package needed to compute and solve system passed on and returned by dynAIMsolver1.
|
||||
( see http://www.federalreserve.gov/Pubs/oss/oss4/aimindex.html )
|
||||
|
||||
The path to the AIM directory is added by dynare_config.m
|
||||
|
||||
DR1 tries to invoke AIM if options_.useAIM == 1 is set and, if not check only, and if 1st order only:
|
||||
if (options_.useAIM == 1) && (task == 0) && (options_.order == 1)
|
||||
|
||||
for start, options_.useAIM = 0 is set by default in global_initialization.m so that system uses mjdgges by default.
|
||||
|
||||
If AIM is to be used, options_.useAIM = 1 needs to be set either in the model <>.mod file, before invoking, estimate and/or stoch_simul, or by issuing appropriate command for estimate and/or stoch_simul.
|
||||
|
||||
NOTE: in the current implementation, as of July 2008, handling of exceptions is rather fundamental and, in particular, when Blanchard and Kahn conditions are not met, only a large penalty value 1.0e+8 is being set.
|
||||
|
||||
Hence, system may not coverge ot resluts may not be accurate if there were many messages like
|
||||
Error in AIM: aimcode=4 : Aim: too few big roots
|
||||
or
|
||||
Error in AIM: aimcode=3 : Aim: too many big roots
|
||||
especially close to the point of convergence.
|
||||
|
||||
However, if other exceptions occur and aimcode (see codes below) is higher than 5, the system resets options_.useAIM = 0 tries to use mjdgges instead.
|
||||
|
||||
|
||||
APPENDIX
|
||||
|
||||
% 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'= £
|
||||
|
||||
|
||||
%function [dr,aimcode]=dynAIMsolver1(jacobia_,M_,dr)
|
||||
% 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.';
|
||||
% 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) for forward looking models, passing into dynAIMsolver
|
||||
% 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
|
||||
% part of DYNARE, copyright Dynare Team (1996-2008)
|
||||
% Gnu Public License.
|
||||
|
|
@ -0,0 +1,20 @@
|
|||
function e = SPAimerr(c);
|
||||
% e = aimerr(c);
|
||||
%
|
||||
% Interpret the return codes generated by the aim routines.
|
||||
|
||||
% The return code c = 2 is used by aim_schur.m but not by aim_eig.m.
|
||||
|
||||
if(c==1) e='Aim: unique solution.';
|
||||
elseif(c==2) e='Aim: roots not correctly computed by real_schur.';
|
||||
elseif(c==3) e='Aim: too many big roots.';
|
||||
elseif(c==35) e='Aim: too many big roots, and q(:,right) is singular.';
|
||||
elseif(c==4) e='Aim: too few big roots.';
|
||||
elseif(c==45) e='Aim: too few big roots, and q(:,right) is singular.';
|
||||
elseif(c==5) e='Aim: q(:,right) is singular.';
|
||||
elseif(c==61) e='Aim: too many exact shiftrights.';
|
||||
elseif(c==62) e='Aim: too many numeric shiftrights.';
|
||||
else e='Aimerr: return code not properly specified';
|
||||
end
|
||||
|
||||
return
|
||||
|
|
@ -0,0 +1,84 @@
|
|||
% [b,rts,ia,nexact,nnumeric,lgroots,aimcode] = ...
|
||||
% SPAmalg(h,neq,nlag,nlead,condn,uprbnd)
|
||||
%
|
||||
% Solve a linear perfect foresight model using the matlab eig
|
||||
% function to find the invariant subspace associated with the big
|
||||
% roots. This procedure will fail if the companion matrix is
|
||||
% defective and does not have a linearly independent set of
|
||||
% eigenvectors associated with the big roots.
|
||||
%
|
||||
% Input arguments:
|
||||
%
|
||||
% h Structural coefficient matrix (neq,neq*(nlag+1+nlead)).
|
||||
% neq Number of equations.
|
||||
% nlag Number of lags.
|
||||
% nlead Number of leads.
|
||||
% condn Zero tolerance used as a condition number test
|
||||
% by numeric_shift and reduced_form.
|
||||
% uprbnd Inclusive upper bound for the modulus of roots
|
||||
% allowed in the reduced form.
|
||||
%
|
||||
% Output arguments:
|
||||
%
|
||||
% b Reduced form coefficient matrix (neq,neq*nlag).
|
||||
% rts Roots returned by eig.
|
||||
% ia Dimension of companion matrix (number of non-trivial
|
||||
% elements in rts).
|
||||
% nexact Number of exact shiftrights.
|
||||
% nnumeric Number of numeric shiftrights.
|
||||
% lgroots Number of roots greater in modulus than uprbnd.
|
||||
% aimcode Return code: see function aimerr.
|
||||
% Modified to deal with infinite or nan A Dec 12 2007
|
||||
% =========================================================
|
||||
function [b,rts,ia,nexact,nnumeric,lgroots,aimcode] = ...
|
||||
SPAmalg(h,neq,nlag,nlead,condn,uprbnd)
|
||||
b=[];rts=[];ia=[];nexact=[];nnumeric=[];lgroots=[];aimcode=[];
|
||||
if(nlag<1 || nlead<1)
|
||||
error('Aim_eig: model must have at least one lag and one lead');
|
||||
end
|
||||
% Initialization.
|
||||
nexact=0;nnumeric=0;lgroots=0;iq=0;aimcode=0;b=0;qrows=neq*nlead;qcols=neq*(nlag+nlead);
|
||||
bcols=neq*nlag;q=zeros(qrows,qcols);rts=zeros(qcols,1);
|
||||
[h,q,iq,nexact]=SPExact_shift(h,q,iq,qrows,qcols,neq);
|
||||
if (iq>qrows)
|
||||
aimcode = 61;
|
||||
return;
|
||||
end
|
||||
[h,q,iq,nnumeric]=SPNumeric_shift(h,q,iq,qrows,qcols,neq,condn);
|
||||
if (iq>qrows)
|
||||
aimcode = 62;
|
||||
return;
|
||||
end
|
||||
[a,ia,js] = SPBuild_a(h,qcols,neq);
|
||||
if (ia ~= 0)
|
||||
if any(any(isnan(a))) || any(any(isinf(a)))
|
||||
disp('A is NAN or INF')
|
||||
aimcode=63;
|
||||
return
|
||||
end
|
||||
[w,rts,lgroots,flag_trouble]=SPEigensystem(a,uprbnd,min(length(js),qrows-iq+1));
|
||||
if flag_trouble==1;
|
||||
disp('Problem in SPEIG');
|
||||
aimcode=63;
|
||||
return
|
||||
end
|
||||
q = SPCopy_w(q,w,js,iq,qrows);
|
||||
end
|
||||
test=nexact+nnumeric+lgroots;
|
||||
if (test > qrows)
|
||||
aimcode = 3;
|
||||
elseif (test < qrows)
|
||||
aimcode = 4;
|
||||
end
|
||||
if(aimcode==0)
|
||||
[nonsing,b]=SPReduced_form(q,qrows,qcols,bcols,neq,condn);
|
||||
if ( nonsing && aimcode==0)
|
||||
aimcode = 1;
|
||||
elseif (~nonsing && aimcode==0)
|
||||
aimcode = 5;
|
||||
elseif (~nonsing && aimcode==3)
|
||||
aimcode = 35;
|
||||
elseif (~nonsing && aimcode==4)
|
||||
aimcode = 45;
|
||||
end
|
||||
end
|
||||
|
|
@ -0,0 +1,39 @@
|
|||
% [a,ia,js] = SPBuild_a(h,qcols,neq)
|
||||
%
|
||||
% Build the companion matrix, deleting inessential lags.
|
||||
% Solve for x_{t+nlead} in terms of x_{t+nlag},...,x_{t+nlead-1}.
|
||||
|
||||
function [a,ia,js] = SPBuild_a(h,qcols,neq)
|
||||
|
||||
left = 1:qcols;
|
||||
right = qcols+1:qcols+neq;
|
||||
hs=SPSparse(h);
|
||||
|
||||
hs(:,left) = -hs(:,right)\hs(:,left);
|
||||
|
||||
% Build the big transition matrix.
|
||||
|
||||
a = zeros(qcols,qcols);
|
||||
|
||||
if(qcols > neq)
|
||||
eyerows = 1:qcols-neq;
|
||||
eyecols = neq+1:qcols;
|
||||
a(eyerows,eyecols) = eye(qcols-neq);
|
||||
end
|
||||
hrows = qcols-neq+1:qcols;
|
||||
a(hrows,:) = hs(:,left);
|
||||
|
||||
% Delete inessential lags and build index array js. js indexes the
|
||||
% columns in the big transition matrix that correspond to the
|
||||
% essential lags in the model. They are the columns of q that will
|
||||
% get the unstable left eigenvectors.
|
||||
|
||||
js = 1:qcols;
|
||||
zerocols = sum(abs(a)) == 0;
|
||||
while( any(zerocols) )
|
||||
a(:,zerocols) = [];
|
||||
a(zerocols,:) = [];
|
||||
js(zerocols) = [];
|
||||
zerocols = sum(abs(a)) == 0;
|
||||
end
|
||||
ia = length(js);
|
||||
|
|
@ -0,0 +1,12 @@
|
|||
% q = SPCopy_w(q,w,js,iq,qrows)
|
||||
%
|
||||
% Copy the eigenvectors corresponding to the largest roots into the
|
||||
% remaining empty rows and columns js of q
|
||||
|
||||
function q = SPCopy_w(q,w,js,iq,qrows)
|
||||
|
||||
if(iq < qrows)
|
||||
lastrows = iq+1:qrows;
|
||||
wrows = 1:length(lastrows);
|
||||
q(lastrows,js) = w(:,wrows)';
|
||||
end
|
||||
|
|
@ -0,0 +1,30 @@
|
|||
function [w,rts,lgroots] = SPEigQZ(a,uprbnd,rowsLeft)
|
||||
% [w,rts,lgroots] = SPEigQZ(a,uprbnd)
|
||||
%
|
||||
% Compute the roots and the left eigenvectors of the companion
|
||||
% matrix, sort the roots from large-to-small, and sort the
|
||||
% eigenvectors conformably. Map the eigenvectors into the real
|
||||
% domain. Count the roots bigger than uprbnd.
|
||||
|
||||
|
||||
eigOpt.disp=0;
|
||||
[w,d] = eig(full(a'));
|
||||
rts = diag(d);
|
||||
mag = abs(rts);
|
||||
[mag,k] = sort(-mag);
|
||||
rts = rts(k);
|
||||
|
||||
ws=SPSparse(w);
|
||||
ws = ws(:,k);
|
||||
|
||||
% Given a complex conjugate pair of vectors W = [w1,w2], there is a
|
||||
% nonsingular matrix D such that W*D = real(W) + imag(W). That is to
|
||||
% say, W and real(W)+imag(W) span the same subspace, which is all
|
||||
% that aim cares about.
|
||||
|
||||
ws = real(ws) + imag(ws);
|
||||
|
||||
lgroots = sum(abs(rts) > uprbnd);
|
||||
|
||||
w=full(ws);
|
||||
|
||||
|
|
@ -0,0 +1,48 @@
|
|||
% [w,rts,lgroots] = SPEigensystem(a,uprbnd)
|
||||
%
|
||||
% Compute the roots and the left eigenvectors of the companion
|
||||
% matrix, sort the roots from large-to-small, and sort the
|
||||
% eigenvectors conformably. Map the eigenvectors into the real
|
||||
% domain. Count the roots bigger than uprbnd.
|
||||
|
||||
function [w,rts,lgroots,flag_trouble] = SPEigensystem(a,uprbnd,rowsLeft)
|
||||
opts.disp=0;
|
||||
try
|
||||
[w,d] = eigs(a',rowsLeft,'LM',opts);
|
||||
rts = diag(d);
|
||||
mag = abs(rts);
|
||||
[mag,k] = sort(-mag);
|
||||
rts = rts(k);
|
||||
catch
|
||||
%disp('Catch in SPE');
|
||||
%pause(0.5);
|
||||
%aStr=datestr(clock);
|
||||
%eval(['save ' regexprep(aStr,' ','') ' a']);
|
||||
try
|
||||
[w,d]=eig(a');
|
||||
catch
|
||||
lasterr
|
||||
w=[];rts=[];lgroots=[];
|
||||
flag_trouble=1;
|
||||
return
|
||||
end
|
||||
rts = diag(d);
|
||||
mag = abs(rts);
|
||||
[mag,k] = sort(-mag);
|
||||
rts = rts(k);
|
||||
end
|
||||
flag_trouble=0;
|
||||
|
||||
ws=SPSparse(w);
|
||||
ws = ws(:,k);
|
||||
|
||||
% Given a complex conjugate pair of vectors W = [w1,w2], there is a
|
||||
% nonsingular matrix D such that W*D = real(W) + imag(W). That is to
|
||||
% say, W and real(W)+imag(W) span the same subspace, which is all
|
||||
% that aim cares about.
|
||||
|
||||
ws = real(ws) + imag(ws);
|
||||
|
||||
lgroots = sum(abs(rts) > uprbnd);
|
||||
|
||||
w=full(ws);
|
||||
|
|
@ -0,0 +1,22 @@
|
|||
% [h,q,iq,nexact] = exact_shift(h,q,iq,qrows,qcols,neq)
|
||||
%
|
||||
% Compute the exact shiftrights and store them in q.
|
||||
|
||||
function [h,q,iq,nexact] = SPExact_shift(h,q,iq,qrows,qcols,neq)
|
||||
|
||||
hs=SPSparse(h);
|
||||
nexact = 0;
|
||||
left = 1:qcols;
|
||||
right = qcols+1:qcols+neq;
|
||||
zerorows = find( sum(abs( hs(:,right)' ))==0 );
|
||||
|
||||
while( any(zerorows) & iq <= qrows )
|
||||
nz = length(zerorows);
|
||||
q(iq+1:iq+nz,:) = hs(zerorows,left);
|
||||
hs(zerorows,:) = SPShiftright(hs(zerorows,:),neq);
|
||||
iq = iq + nz;
|
||||
nexact = nexact + nz;
|
||||
zerorows = find( sum(abs( hs(:,right)' ))==0 );
|
||||
end
|
||||
h=full(hs);
|
||||
|
||||
|
|
@ -0,0 +1,26 @@
|
|||
% [h,q,iq,nnumeric] = ...
|
||||
% SPNumeric_shift(h,q,iq,qrows,qcols,neq,condn)
|
||||
%
|
||||
% Compute the numeric shiftrights and store them in q.
|
||||
|
||||
function [h,q,iq,nnumeric] = SPNumeric_shift(h,q,iq,qrows,qcols,neq,condn)
|
||||
|
||||
nnumeric = 0;
|
||||
left = 1:qcols;
|
||||
right = qcols+1:qcols+neq;
|
||||
|
||||
[Q,R,E] = qr( h(:,right) );
|
||||
zerorows = find( abs(diag(R)) <= condn );
|
||||
|
||||
while( any(zerorows) & iq <= qrows )
|
||||
h=sparse(h);
|
||||
Q=sparse(Q);
|
||||
h = Q'*h;
|
||||
nz = length(zerorows);
|
||||
q(iq+1:iq+nz,:) = h(zerorows,left);
|
||||
h(zerorows,:) = SPShiftright( h(zerorows,:), neq );
|
||||
iq = iq + nz;
|
||||
nnumeric = nnumeric + nz;
|
||||
[Q,R,E] = qr( full(h(:,right)) );
|
||||
zerorows = find( abs(diag(R)) <= condn );
|
||||
end
|
||||
|
|
@ -0,0 +1,47 @@
|
|||
% scof = SPObstruct(cof,cofb,neq,nlag,nlead)
|
||||
%
|
||||
% Construct the coefficients in the observable structure.
|
||||
%
|
||||
% Input arguments:
|
||||
%
|
||||
% cof structural coefficients
|
||||
% cofb reduced form
|
||||
% neq number of equations
|
||||
% nlag number of lags
|
||||
% nlead number of leads
|
||||
%
|
||||
% Output arguments:
|
||||
%
|
||||
% scof observable structure coefficients
|
||||
|
||||
function scof = SPObstruct(cof,cofb,neq,nlag,nlead)
|
||||
|
||||
|
||||
% Append the negative identity to cofb
|
||||
|
||||
cofb = [cofb, -eye(neq)];
|
||||
scof = zeros(neq,neq*(nlag+1));
|
||||
q = zeros(neq*nlead, neq*(nlag+nlead));
|
||||
[rc,cc] = size(cofb);
|
||||
qs=sparse(q);
|
||||
qs(1:rc,1:cc) = sparse(cofb);
|
||||
qcols = neq*(nlag+nlead);
|
||||
|
||||
if( nlead > 1 )
|
||||
for i = 1:nlead-1
|
||||
rows = i*neq + (1:neq);
|
||||
qs(rows,:) = SPShiftright( qs((rows-neq),:), neq );
|
||||
end
|
||||
end
|
||||
|
||||
l = (1: neq*nlag);
|
||||
r = (neq*nlag+1: neq*(nlag+nlead));
|
||||
|
||||
qs(:,l) = - qs(:,r) \ qs(:,l);
|
||||
|
||||
minus = 1: neq*(nlag+1);
|
||||
plus = neq*(nlag+1)+1: neq*(nlag+1+nlead);
|
||||
|
||||
cofs=sparse(cof);
|
||||
scof(:,neq+1:neq*(nlag+1)) = cofs(:,plus)*qs(:,l);
|
||||
scof = scof + cofs(:,minus);
|
||||
|
|
@ -0,0 +1,91 @@
|
|||
function [q,rts,ia,nexact,nnumeric,lgroots,aimcode] = ...
|
||||
SPQZ(aa,bb)
|
||||
% [q,rts,ia,nexact,nnumeric,lgroots,aimcode] = ...
|
||||
% SPQZ(aa,bb)
|
||||
%
|
||||
% Solve a linear perfect foresight model using the matlab eig
|
||||
% function to find the invariant subspace associated with the big
|
||||
% roots. This procedure will fail if the companion matrix is
|
||||
% defective and does not have a linearly independent set of
|
||||
% eigenvectors associated with the big roots.
|
||||
%
|
||||
% Input arguments:
|
||||
%
|
||||
% h Structural coefficient matrix (neq,neq*(nlag+1+nlead)).
|
||||
% neq Number of equations.
|
||||
% nlag Number of lags.
|
||||
% nlead Number of leads.
|
||||
% condn Zero tolerance used as a condition number test
|
||||
% by numeric_shift and reduced_form.
|
||||
% uprbnd Inclusive upper bound for the modulus of roots
|
||||
% allowed in the reduced form.
|
||||
%
|
||||
% Output arguments:
|
||||
%
|
||||
% b Reduced form coefficient matrix (neq,neq*nlag).
|
||||
% rts Roots returned by eig.
|
||||
% ia Dimension of companion matrix (number of non-trivial
|
||||
% elements in rts).
|
||||
% nexact Number of exact shiftrights.
|
||||
% nnumeric Number of numeric shiftrights.
|
||||
% lgroots Number of roots greater in modulus than uprbnd.
|
||||
% aimcode Return code: see function aimerr.
|
||||
[neq,aCols]=size(aa);
|
||||
[bRows,bCols]=size(bb);
|
||||
if(~and(neq ==aCols,and(neq== bRows,neq ==bCols)))
|
||||
error(SPQZ:'aa, bb must be square and same dim')
|
||||
end
|
||||
nlag=0
|
||||
nlead=1
|
||||
b=[];
|
||||
rts=[];
|
||||
ia=[];
|
||||
nexact=[];
|
||||
nnumeric=[];
|
||||
lgroots=[];
|
||||
aimcode=[];
|
||||
|
||||
condn=1.0e-8
|
||||
uprbnd=1.0e16
|
||||
|
||||
% Initialization.
|
||||
nexact = 0;
|
||||
nnumeric = 0;
|
||||
lgroots = 0;
|
||||
iq = 0;
|
||||
aimcode = 0;
|
||||
b=0;
|
||||
qrows = neq*nlead;
|
||||
qcols = neq*(nlag+nlead);
|
||||
bcols = neq*nlag;
|
||||
q = zeros(qrows,qcols);
|
||||
rts = zeros(qcols,1);
|
||||
|
||||
% Compute the auxiliary initial conditions and store them in q.
|
||||
h=[aa bb];
|
||||
|
||||
|
||||
[h,q,iq,nexact] = SPExact_shift(h,q,iq,qrows,qcols,neq);
|
||||
if (iq>qrows)
|
||||
aimcode = 61;
|
||||
return;
|
||||
end
|
||||
|
||||
[h,q,iq,nnumeric] = SPNumeric_shift(h,q,iq,qrows,qcols,neq,condn);
|
||||
if (iq>qrows)
|
||||
aimcode = 62;
|
||||
return;
|
||||
end
|
||||
|
||||
% Build the companion matrix. Compute the stability conditions, and
|
||||
% combine them with the auxiliary initial conditions in q.
|
||||
|
||||
[a,ia,js] = SPBuild_a(h,qcols,neq);
|
||||
|
||||
if (ia ~= 0)
|
||||
[w,rts,lgroots] = SPEigQZ(a,uprbnd,min(length(ia),qrows-iq+1));
|
||||
%enough in w to fill q
|
||||
q = SPCopy_w(q,w,js,iq,qrows);
|
||||
end
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,26 @@
|
|||
% [nonsing,b] = SPReduced_form(q,qrows,qcols,bcols,neq,b,condn);
|
||||
%
|
||||
% Compute reduced-form coefficient matrix, b.
|
||||
|
||||
function [nonsing,b] = SPReduced_form(q,qrows,qcols,bcols,neq,condn);
|
||||
b=[];
|
||||
qs=SPSparse(q);
|
||||
left = 1:qcols-qrows;
|
||||
right = qcols-qrows+1:qcols;
|
||||
nonsing = rcond(full(qs(:,right))) > condn;
|
||||
if(nonsing)
|
||||
qs(:,left) = -qs(:,right)\qs(:,left);
|
||||
b = qs(1:neq,1:bcols);
|
||||
b = full(b);
|
||||
else %rescale by dividing row by maximal qr element
|
||||
%'inverse condition number small, rescaling '
|
||||
themax=max(abs(qs(:,right)),[],2);
|
||||
oneover = diag(1 ./ themax);
|
||||
nonsing = rcond(full(oneover *qs(:,right))) > condn;
|
||||
if(nonsing)
|
||||
qs(:,left) = -(oneover*qs(:,right))\(oneover*(qs(:,left)));
|
||||
b = qs(1:neq,1:bcols);
|
||||
b = full(b);
|
||||
end
|
||||
end
|
||||
|
||||
|
|
@ -0,0 +1,16 @@
|
|||
function [y] = SPShiftright(x,n)
|
||||
|
||||
% [y] = shiftright(x,n)
|
||||
%
|
||||
% Shift the rows of x to the right by n columns, leaving zeros in the
|
||||
% first n columns.
|
||||
|
||||
[rows,cols] = size(x);
|
||||
|
||||
left = 1:cols-n;
|
||||
right = n+1:cols;
|
||||
|
||||
y = zeros(rows,cols);
|
||||
y(:,right) = x(:,left);
|
||||
|
||||
return
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
function res = SPSparse(aMat)
|
||||
res=sparse(aMat);
|
||||
|
|
@ -0,0 +1,118 @@
|
|||
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
|
||||
% 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.';
|
||||
% 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
|
||||
% part of DYNARE, copyright Dynare Team (1996-2008)
|
||||
% Gnu Public License.
|
||||
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(:)));
|
||||
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,ia,nexact,nnumeric,lgroots,aimcode] =...
|
||||
SPAmalg(theAIM_H,neq, lags,leads,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
|
||||
dr.ghx= bb_ord(:,find(any(bb_ord,1))); % get non-zero columns only
|
||||
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);
|
||||
scof1= scof1(:,dr.order_var);
|
||||
dr.ghu =scof1\theAIM_Psi;
|
||||
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
|
||||
Loading…
Reference in New Issue