@--This program good for bond economy

Decide what cutoff value to use to distinguish unstable from stable eigenvalues--@

cutoff=1.00001;
critc=0.0000001;


dlibrary -a PszTgSen;
#include PsTgSen.prc;

lhs=(gradp(&lhs1,xbar)|gradp(&lhs2,xbar)|gradp(&lhs3,xbar)|
gradp(&lhs4,xbar)|gradp(&lhs5,xbar)|gradp(&lhs6,xbar)|
gradp(&lhs7,xbar)|gradp(&lhs8,xbar)|gradp(&lhs9,xbar)|
gradp(&lhs10,xbar)|gradp(&lhs11,xbar)|gradp(&lhs12,xbar));

rhs=(gradp(&rhs1,xbar)|gradp(&rhs2,xbar)|gradp(&rhs3,xbar)|
gradp(&rhs4,xbar)|gradp(&rhs5,xbar)|gradp(&rhs6,xbar)|
gradp(&rhs7,xbar)|gradp(&rhs8,xbar)|gradp(&rhs9,xbar)|
gradp(&rhs10,xbar)|gradp(&rhs11,xbar)|gradp(&rhs12,xbar));


@--Now use Paul Soderlind's procedure to give alternative representations for lhs and rhs.
Use Golub and Van Loan notation: lhs=B=QSZh, rhs=A=QTZh--@

{ t,s,q,z, lamdap } = zTgSen( rhs,lhs,cutoff,1,0);

ns=rows(selif(abs(lamdap),abs(lamdap).<cutoff));

if ns>=nvar;

@--All eigenvalues are less than cutoff - no unstable roots - nothing more to do--@

elseif ns<nvar;

@--Partition the Q, Z, S and T matrices--@

qh=inv(q);

qh11=qh[1:ns,1:ns];
qh12=qh[1:ns,ns+1:nvar];
qh21=qh[ns+1:nvar,1:ns];
qh22=qh[ns+1:nvar,ns+1:nvar];

z11=z[1:ns,1:ns];
z12=z[1:ns,ns+1:nvar];
z21=z[ns+1:nvar,1:ns];
z22=z[ns+1:nvar,ns+1:nvar];

s11=s[1:ns,1:ns];
s12=s[1:ns,ns+1:nvar];
s21=s[ns+1:nvar,1:ns];
s22=s[ns+1:nvar,ns+1:nvar];

t11=t[1:ns,1:ns];
t12=t[1:ns,ns+1:nvar];
t21=t[ns+1:nvar,1:ns];
t22=t[ns+1:nvar,ns+1:nvar];

@--The following terms are defined as follows.

x(t)[ns+1:nvar] = FA * x(t)[1:ns]

f(t)[ns+1:nvar] = FC * f(t)[1:ns]

x(t+1)[1:ns] = FB * x(t)[1:ns] + FD * f(t+1)[1:ns]

f(t+1)= 0 | 0 | 0 | shock1 | shock2 --@

FA=z21*inv(z11);
FB=z11*inv(s11)*t11*inv(z11);
FC=-inv(qh22)*qh21;
FD=z11*inv(s11)*(qh11+qh12*FC);

print "EIGENVALUES";
print "lamdap " lamdap;
print "###################";

print "DECISION RULES";
print "fa " fa;
print "fb " fb;
print "fd " fd;
print "###################";


endif;

dlibrary -d;