function [M,N] = ComAlgS( A,B,Q,R,U,bet,cutoff );
%ComAlgS    Solves for optimization problem in commitment case, stable system.
%
%
%
%  Usage:     [M,N] = ComItAlg( A,B,Q,R,U,bet,cutoff );
%
%  Input:     A        n1xn1 matrix
%             B        n1xk matrix
%             Q        n1xn1 matrix, symmetric
%             R        kxk matrix, symmetric
%             U        n1xk matrix
%             bet      scalar, discount factor (eg 0.99)
%             cutoff   scalar, abs(roots) > cutoff indicates instability
%
%  Output:    M        n1xn1 matrix, x1(t+1) = M*x1(t) + e(t+1)
%             N        (k+n1)xn1 matrix, [u(t),p1(t)] = N*x1(t)
%
%
%
%  Details:   Solve for optimization problem in commitment case. The economy
%             evolves as
%
%             x1(t+1)       =   A * x1(t)    + Bu(t) + e(t+1)
%
%             where x1(t) ("backward looking") has n1 elements. The initial
%             value of x1, x1(0), is given by history.
%             u(t) is a kx1 vector of decision variables, which follows
%             the decision rule.
%
%             The loss function of the policy maker is
%
%
%             Sum{ (bet^t)*[x1(t)'Q*x1(t)+2x1(t)'U*u(t)+u(t)'R*u(t)],t=0,1,... }
%
%
%  The solution is summarized by the matrices M and N in
%
%  x1(t+1)       = M*x1(t) + e(t+1)
%  [u(t),p1(t)]  = N*x1(t)
%
%
%  Calls on:   reorder (if MatLab)
%
%
%  Paul Söderlind, Paul.Soderlind@unisg.ch, June 2001
%----------------------------------------------------------------------------

Q = (Q + Q')/2;                %to make symmetric
R = (R + R')/2;

n1 = size(A,1);                %no of variables (all stable)
k  = size(R,1);                %no. of control variables

G =  [ eye(n1),       zeros(n1,k),   zeros(n1,n1) ;
       zeros(n1,n1),  zeros(n1,k),   (bet*A')     ;
       zeros(k,n1),   zeros(k,k),    (-B')        ];

D =  [ A,             B,             zeros(n1,n1) ;
       (-bet*Q),     (-bet*U),       eye(n1)      ;
       U',            R,             zeros(k,n1)  ];


if exist('OCTAVE_VERSION');         %Octave, real generalized Schur, odd def of lambda
  [S,T,Z,lambda] = qz(G,D,'B');
else;                                    %MatLab
  [S,T,Qa,Z] = qz(G,D);    %MatLab: G=Q'SZ' and D=Q'TZ'; Paul S:  G=QSZ' and D=QTZ', but Q isn't used
  [S,T,Qa,Z] = reorder(S,T,Qa,Z);   % reordering of generalized eigenvalues in ascending order
end;


Stt = S(1:n1,1:n1);
Ttt = T(1:n1,1:n1);
Zkt = Z(1:n1,1:n1);
Zlt = Z(n1+1:n1+k+n1,1:n1);

if cond(Zkt) > 1e+14;
  warning('Zkt is singular: rank condition for solution not satisfied');
  M = NaN; N = NaN;
  return;
end;

Zkt_1 = inv(Zkt);         %inverting
Stt_1 = inv(Stt);


M = real(Zkt*Stt_1*Ttt*Zkt_1);      %x1(t+1)      = M*x1(t)+e(t+1)
N = real(Zlt*Zkt_1);                %[u(t),p1(t)] = N*x1(t)
                                    %getting rid of any trivial complex parts

maxEig = max(abs(eig(M)));
if maxEig > cutoff;
  warning('Some abs(eigenvalue) of M > 1');
end;
%-----------------------------------------------------------------------