function [M,N,Ma,Na,Fv,J0,J0a] = ComAlgSR( A,B,Q,R,U,bet,cutoff,C,theta );
%ComAlgSR    Solves for robust optimization problem in commitment case, for a
%            backward-looking system (no forward-looking variables).
%
%            The economy evolves according to
%
%            x(t+1) = A * x(t) + Bu(t) + C*e(t+1),
%
%            where x(t) is nx1, u(t) is kx1 and e(t+1) is nx1
%
%            The loss function of the policy
%            maker is
%
%            Sum{ (bet^t)*[x(t)'Q*x(t)+2x(t)'U*u(t)+u(t)'R*u(t)],t=0,1,... }
%
%
%
%
%  Usage:    [M,N,Ma,Na,Fv,J0,J0a] = ComAlgSR( A,B,Q,R,U,bet,cutoff,C,theta );
%
%
%  Input:    A        see ComAlgR
%            B
%            Q
%            R
%            U
%            bet
%            cutoff
%            C
%            theta
%
%  Output:   M        nxn matrix, x(t+1) = M*x(t) + C*e(t+1), worst case
%            N        (k+n*2)x(n) matrix, [u;v;p1]= N*x, worst case
%            Ma       nxn matrix, x(t+1) = Ma*x(t) + C*e(t+1), approximating model
%            Na       (k+n*2)x(n) matrix, [u;v;p1]= Na*x, approximating model
%            Fv       nxn matrix, v = -Fv*x
%            J0       loss function value, worst case
%            J0a      loss function value, approximating model
%
%
%  Note:    N and Na differs only in that rows corrsponding to v are zero in Na.
%
%
%  Calls on: ComAlgS
%
%  Paolo Giordani and Paul.Soderlind@unisg.ch, Jun 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


Btilde = [ B, C ];
Utilde = [ U, zeros(n1,n1) ];
Rtilde = [ R,            zeros(k,n1); ...
           zeros(n1,k),  -theta*eye(n1) ];


[M,N] = ComAlgS( A,Btilde,Q,Rtilde,Utilde,bet,cutoff );   %commitment, robust policy


N2 = N(1:k,:);                        % u = N2*x1; note: the P matrix is eye(n1) since n2=0
N3 = N(k+1:k+n1,:);                   % v = N3*x1

Fu  = -N2;                       %u = -Fu*x1
Fv  = -N3;                       %v = -Fv*x1

Mw = A-B*Fu-C*Fv;     %alternative expression for M, should equal M
%disp(' ');disp('M-Mw');disp([M-Mw]);


Ma = A-B*Fu;                     %x1(t+1) = Ma*x1(t), approximating model
Na = N;                          %[u;v;p1] = Na*x1, approximating model
Na(k+1:k+n1,:) = zeros(n1,n1);   %no evil shocks in approximate model

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



Ptilde = [ eye(n1); -Fu; -Fv ];             %value function, worst case
W      = Ptilde' * [Q, Utilde; Utilde', Rtilde] * Ptilde;
vecV   = (eye(n1*n1) - kron(M',bet*M'))\W(:);
V      = reshape(vecV,n1,n1);                 %solves V = W + bet*M'*V*M
if bet == 1;
  J0 = trace(C'*V*C*eye(n1));
else;
  J0 = bet/(1-bet) * trace(C'*V*C*eye(n1));
end;


Ptilde = [ eye(n1); -Fu; zeros(n1,n1) ];             %value function
W      = Ptilde' * [Q, Utilde; Utilde', Rtilde] * Ptilde;
vecV   = (eye(n1*n1) - kron(Ma',bet*Ma'))\W(:);
V      = reshape(vecV,n1,n1);                 %solves V = W + bet*Ma'*V*Ma
if bet == 1;
  J0a = trace(C'*V*C*eye(n1));
else;
  J0a = bet/(1-bet) * trace(C'*V*C*eye(n1));
end;
%-----------------------------------------------------------------------


% OLD STUFF

% Vs_1   = 0;
% Vb     = 0;
% Vdiff  = 100;
% while Vdiff > 1E-6;    %iterating on value matrix V
%   Vb    = W + bet*M'*Vs_1*M;  %J(0) = Sum[(bet^t)z(t)'Wz(t),t=0,...]
%   Vdiff = max(  max( abs(Vs_1-Vb) )  );
%   Vs_1  = Vb;            %Vs_1 is V(s-1), V is V(s)
% end;
% disp('testing V ');disp([V,V-Vb]);

% Vs_1   = 0;
% Vb     = 0;
% Vdiff  = 100;
% while Vdiff > 1E-6;
%   Vb    = W + bet*Ma'*Vs_1*Ma;
%   Vdiff = max(  max( abs(Vs_1-Vb) )  );
%   Vs_1  = Vb;
% end;
% disp('testing V ');disp([V,V-Vb]);