function OUT = sar(y,X,W,eigvals)

%Version of SAR for simulation output


%Written by: TONY E. SMITH, 4/10/98 (Modified: 4/17/00)

%NOTE: P-values are only reported if StatToolbox is running.


%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: OUT = nX4 matrix with columns:
%        OUT(:,1) = SAR predicted yHats .
%        OUT(:,2) = Observed y values .
%        OUT(:,3) = SAR residual vector, r .
%        OUT(:,4) = Weighted residuals, W*r .
%        OUT(:,5) = OLS residuals, r0 .
%        OUT(:,6) = Weighted OLS residuals, W*r0 .
%
%               I = inverse of Asymptotic covariance matrix
%               
%        NOTE: The residual vector, r, is for the reduced
%              OLS model, B*y = B*X*b + e. Hence, a comparison
%              of the regressions, (r0 on W*r0) and (r on W*r),
%              indicates the degree to which the SAR model has
%              reduced (W-type) spatial autocorrelation.
%
%
%SCREEN OUTPUT:
%
%         TABLE1: (REGRESSION RESULTS)
%                 beta coeffs plus z-values for significance tests
%                 (If MATLAB Statistics Toolbox is present, the P-values
%                  are also displayed as 'PROB').
%
%         TABLE2: (GOODNESS OF FIT)
%                 Ps_R2 = Pseudo R-Square (in reduced model)
%                 Sq-Corr = squared correlation coefficient
%                 Lik = value of log likelihood
%
%         TABLE3: (AUTOCORRELATION PARAMETER)
%                 Rho estimate plus z-value for significance tests
%                 (If MATLAB Statistics Toolbox is present, the P-values
%                  are also displayed as 'PROB').
%
%         TABLE4: (TESTS OF SPATIAL AUTOCORRELATION MODEL)
%                 LR  = likelihood ratio for spatial autocorrelation
%                 Com-LR = likelihood ratio for common-factor hyp.

%PLOT OUTPUT: See the program SAR_PLOT to plot output.



[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 ; %OLS residuals

var0 = (1/n)*r0'*r0 ; %OLS error variance

%Eigen values

if nargin == 5
   
   val = real(eigvals);
   
   v0 = min(val);
   
   v1 = max(val);
   
elseif nargin == 4
   
   val = real(eig(full(W)));

	[II,JJ,V] = find(val); %find nonzero values

	v0 = min(V);

	v1 = max(V) ;

else
   
   error('Wrong number of input arguments');
   
end


low = max(-1,1/(-abs(v0))) ;

high = 1/v1 ;


%Estimate rho and beta

criterion = .0001; %internal convergence critirion
itermax = 100;     %internal iteration limit
converge = 1;

econverge = r0;
iter = 1;
e = r0;

while (converge > criterion & iter < itermax) 
   
   rho = fmin('sar_lik_1',low,high,[],low,high,val,W,e,X) ;
   
	B = speye(n) - rho * sparse(W) ;

	b = inv(X'*B'*B*X)*X'*B'*B*y;

	e = y - X * b ;  
   
   converge = sum(abs(e - econverge));

	econverge = e ;

	iter = iter + 1;

end %end estimation loop

if (iter == itermax)
   
   error('No convergence within 100 iterations');
   
end

%Compute Log Likelihood 

var = (1/n) * e'*B'*B*e ; %Final residual variance

sig = sqrt(var);

u = ones(length(val),1) ;

LogDet = u'*log(1 - rho*val) ;

L = -(n/2)*(1 + log(2*pi) + log(var)) + LogDet ;

LR = n * (log(var0) - log(var)) + 2 * LogDet ;

%Prepare output

OUT = [b',rho,sig];