solve_stochastic(A,B,V,W,R,method);

function [S,Q,N1,N2,info] = solve_stochastic(A,B,V,W,R,method);
%  PURPOSE: solves stochastic part of the mdel 
%
%               0 = A y_t + B y_t+1 + C E_t{y_t+1} + V epsilon_t + W epsilon_t+1
%
%           where epsilon_t+1 is i.i.d. random variable with zero mean
%           and holds the transversality condition 
%
%               lim (t->infty) xi^t E_0 y_t = 0
%
%           the model is represented in the state-space form
%
%               y_t     = R u_t   + S_1 epsilon_t + S_2 omega_t
%               u_t     = P u_t-1 + Q_2 epsilon_t + Q_2 omega_t
%
%           where u_t is a state variable, omega_t is i.i.d. random
%           variable with mean zero, X1, X2, Y1, Y2 are any
%           matrices with appropriate size.
%
%           Matrices S, N1, Q, N2 satisfy
%
%               A S_1 + V  = 0          B (R Q_1 + S_1) + W = 0
%               A S_2      = 0          B (R Q_2 + S_2)     = 0
%
%           if matrices S_1, S_2, Q_1, Q_2 exist, then 
%
%               S_1 = S + N1 X_1        Q_1 = Q + N2 X_1
%               S_2 =     N1 X_2        Q_2 =     N2 X_2
%
%           for any matrices X_1, X_2
%           
% ---------------------------------------------------
%  USAGE: [S,Q,N1,N2,info] = solve_stochastic(A,B,V,W,R,method);
%  where: 
%         A,B,V,W                   matrices representing the model,
%         R                         a matrix representing solution to the
%                                   model
%         method                    method to calculate null spaces
%                                   (optional)
%                                       1 - qr with pivoting (default)
%                                       2 - svd
%
%         S,Q,N1,N2                 matrices representing solution to the
%                                   stochastic part of the model
%         info                      return 1 is solution exists and zero
%                                   otherwise. If number of output
%                                   variables is lower than five, then
%                                   an error is thrown if there are no
%                                   solutions.
%
%   COMMENTS:
%
% Copyright  (c) Pawel Kowal (2006)
% All rights reserved
% LREM_SOLVE toolbox is available free for noncommercial academic use only.
% [e-mail]

if nargin<6
    method      = 1;
end
if nargout==5
    throw_error = false;
else
    throw_error = true;
end
info            = 1;
tol             = 100*max(size(A)')*norm(A,1)*eps;

[M1,M2,info1]       = lineareq(A,V,method,tol);
if info1==0
    if throw_error
        error('there are no solutions');
    else
        S       =[];
        Q       =[];
        N1      =[];
        N2      =[];
        info    = 0;
        return;
    end      
end

[M3,M4,info2]       = lineareq(B*[R M2],B*M1+W,method,tol);
if info2==0
    if throw_error
        error('there are no solutions');
    else
        S       =[];
        Q       =[];
        N1      =[];
        N2      =[];
        info    = 0;
        return;
    end      
end

n                   = size(R,2);
Q                   = M3(1:n,:);
N2                  = M4(1:n,:);
S                   = M1 + M2*M3(n+1:end,:);
N1                  = M2*M4(n+1:end,:);