/*
** Written by Paul Klein. Last modified 7 April 1998. 
**
** Format: {ff,pp} = solve(aa,bb,cc,dd,phi,nk);
**
** The meaning of these symbols is explained below.
**
** Suppose you have the system
**
** aa*E[x(t+1)|F(t)]+bb*x(t)+cc*E[z(t+1)|F(t)]+dd*z(t) = 0
**
** where z has the dynamics
**
** z(t+1)=phi*z(t)+eps(t+1)
**
** Partition x(t)=(d(t)'~k(t)')' and let k(t) be predetermined (have zero
** one-period-ahead expectations error).
**
** We assume that eps is an exogenous white noise process and F is the 
** filtration generated by eps, k0 and z0.
**
** Note that k0 and z0 are exogenously given, but that d0 is determined
** endogenously (so as to ensure stability).
**
** Note that nk is the dimension of k(t).
**
** Let w(t) = (k(t)'~z(t)')'
**
** Then this function delivers the matrices ff and pp in the following 
** representation of the unique stable solution (supposing there is one).
**
** d(t)   = ff*w(t)
**
** w(t+1) = pp*w(t) + (0'~eps[t+1]')'
*/

proc (2) = solve(aa,bb,cc,dd,phi,nk);
local nx,nz,nw,ff,pp,s,t,z,q,nd;
local s22,t22,z22,z12,z22i,s22i;
local flag,lambda;

nx = rows(aa);
nz = rows(phi);

dd=-dd;
bb=-bb;

aa = aa~cc|zeros(nz,nx)~eye(nz);
bb = bb~dd|zeros(nz,nx)~phi;

nw = nk+nz;
nx = nx+nz;
nd = nx-nw;

/* Doing the sorted generalized Schur decomposition */

{s,t,z,q,lambda} = psztgsen(aa',bb',1,1,0);

s = s';
t = t';

/* At this stage we have q*a = s*z', q*b = t*z' with q,z unitary and s,t lower triangular. Mean-
   while, the eigenvalue pairs s[i,i],t[i,i] are ordered so that those with 
   abs(s[i,i])<abs(t[i,i]) (the unstable ones) are in the upper left-hand block, and those with
   abs(s[i,i])>abs(t[i,i]) (the stable ones) are in the lower right-hand block.  */

lambda = diag(s)~diag(t);
flag = varput(lambda,"lambda");

/* Checking whether the no of unstable eigenvalues is equal to nd  */

if (abs(s[nd,nd])>abs(t[nd,nd])) or (abs(s[nd+1,nd+1])<abs(t[nd+1,nd+1]));
   print "Warning: Wrong number of stable eigenvalues.";
endif;

z12  = z[1:nd,nd+1:nx];
z22  = z[nd+1:nx,nd+1:nx];

/* Checking whether zh22 is invertible */

if rank(z22)<rows(z22);
print "Fatal error: Invertibility condition violated.";
end;
endif;

z22i = eye(nw)/z22;
s22  = s[nd+1:nx,nd+1:nx];
s22i = eye(nw)/s22;
t22  = t[nd+1:nx,nd+1:nx];

ff  = z12*z22i;
pp  = z22*s22i*t22*z22i;

retp(ff,pp);
endp;