function [b,lam,var] = sp_lag_silent(y,X,W,E)

%This version of SP_LAG runs silently with no screen outputs, 
%and yields only bhat and lambda.


%SP_LAG computes the Spatial Lag Model of y on X given lag weights
%    matrix,W (and using the multivariate FMINS procedure in MATLAB).

%Written by: TONY E. SMITH, 4/11/98 (Modified 6/14/99)

%INPUTS: (i)   y = n-vector of Yi values
%        (ii)  X = (nxk)-matrix of explanatory variables
%                  (Xi1,..,Xik) [NOT including unit vector]
%        (iii) W = (nxn)-matrix of spatial weights
%        (iv)  vnames = names of variables in X. If the variables 
%              are say "altitude" and "temperature", then set:
%              vnames = strvcat('altitude','temperature');For the 
%              case of a single variable, say "altitude", just 
%              write 'altitude'.
%        (v)   eigvals (optional) = vector if eigenvalues of W.
%              [NOTE: for large W-matrices it pays to compute
%               eigenvalues ONCE and then use these as inputs]
%
%OUTPUT: b = SP_LAG regression coefficients


[n,k] = size(X) ;

%Augment X with unit vector

% X0 = X ;
% 
% clear('X') ;
% 
% X(:,1) = ones(n,1) ;
% 
% X(:,2:k+1) = X0 ;

%Initial OLS estimation of beta

b = X\y ;

r0 = y - X * b ;

var0 = (1/n)*r0'*r0 ; %OLS error variance

%Eigen values

if nargin == 4
    
    val = E;
    
    v0 = min(E);

    v1 = max(E) ;
    
    if v0 ~= 0
        
        low = 1/(-abs(v0)) ;
        
    else 
        
        low = 0;
        
    end

    high = 1/v1 ;
    
else
    
    val = real(eig(full(W)));

    [II,JJ,V] = find(val); %find nonzero values

    v0 = min(V);

    v1 = max(V) ;

    if v0 ~= 0
        
        low = 1/(-abs(v0)) ;
        
    else 
        
        low = 0;
        
    end

    high = 1/v1 ;

end

%Initial Estimate of lamda

M = eye(n) - X*inv(X'*X)*X' ;

e0 = M*y ;

e1 = M*(W*y) ;


z = fminbnd('sp_lag_lik',low,high,[],low,high,val,e0,e1,n) ; 


%Compute Log Likelihood

lam = min([z ,high - .001]) ;

lam = max([lam,low + .001]) ;

W = sparse(W);

B = speye(n) - lam*W;

b = inv(X'*X)*X'*B*y ;

e = y - lam*W*y - X*b ;

var = (1/n) * e'*e ;




