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

%SAR computes Simultaneous AutoRegression of y on X given weights W 
%    using the multivariate FMINS procedure in MATLAB.(NOTE: This 
%    program assumes that Statistical Toolbox is running.)

%Written by: TONY E. SMITH, 4/10/98

%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)  varargin = names of variables in X. If the variables 
%              are say "altitude" and "temperature", then the inputs
%              are: (y,X,L,'altitude','temperature')
%
%OUTPUT: OUT = nX4 matrix with columns:
%        OUT(:,1) = OLS residual vector, r0 .
%        OUT(:,2) = Weighted residuals, W*r0 .
%        OUT(:,3) = SAR residual vector, r .
%        OUT(:,4) = Weighted residuals, W*r .
%
%        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.


[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

%Compute eigen values

val = eig(full(W));

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

v0 = min(V);

v1 = max(V) ;

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

high = 1/v1 ;


%Estimate rho
 

rho = fmin('sar_lik',low,high,[],low,high,val,W,y,X) ; 
   

%Compute Log Likelihood

B = sparse(eye(n) - rho * W) ;

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

e = y - X * b ;  

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

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 ;


%COMMON FACTOR HYPOTHESIS

XX(:,1:k+1) = X ;

XX(:,(k+2):(2*k)+1) = W*X0 ; %Note that X0 is used to avoid (1,W*1) singularity

if rank(XX'*XX) == (2*k)+1
         
   M = eye(n) - XX*inv(XX'*XX)*XX' ;
   
   e0 = M*y ;
   
   e1 = M*(W*y) ;
   
   z = fmin('sp_lag_lik',low,high,[],low,high,val,e0,e1) ; 
   
   
   %Compute Log Likelihood
   
   lam = min([z ,high - .001]) ;
   
   lam = max([lam,low + .001]) ;
   
   B0 = eye(n) - lam * W;
   
   b0 = inv(XX'*XX)*XX'*B0*y ;
   
   e0 = y - lam*W*y - XX*b0 ;
   
   var0 = (1/n) * e0'*e0 ;
   
   u = ones(length(val),1);
   
   LogDet = u'*log(1 - lam*val) ;
   
   L0 = -(n/2)*(1 + log(2*pi) + log(var0)) + LogDet ;
   
   LR0 = 2*(L0 - L) ;
   
else
   
   Comment = {'Common Factor Test involves a SINGULARITY'}
   
end

%Calculate Asymptotic Covariance Matrix for the parameter
%vector, (b,var,rho)

G = W*inv(sparse(B)) ;

I(1:k+1,1:k+1) = (1/var)*X'*B'*B*X ;

I(1:k+1,k+2) = zeros(k+1,1) ;

I(1:k+1,k+3) = zeros(k+1,1) ;

I(k+2,1:k+1) = I(1:k+1,k+2)' ;

I(k+2,k+2) = n/(2*(var^2)) ;

I(k+2,k+3) = (1/var)*trace(G) ;

I(k+3,1:k+1) = I(1:k+1,k+3)' ;

I(k+3,k+2) = I(k+2,k+3) ;

I(k+3,k+3) = trace(G*(G + G')) ;

cov = inv(I) ;

se = sqrt(diag(cov)) ;


%SCREEN OUTPUT

clear('v','C');

%Add 'const'

v(1,1) = {'const'};

v(1,2:1+length(varargin)) = varargin;

%Make TABLE1

C(:,1) = v ;

C(:,2) = num2cell(b);

z = abs(b./se(1:k+1)) ;

C(:,3) = num2cell(z);

if exist('normcdf') == 2  %check to see if statistics toolbox exists
   
   C(:,4) = num2cell(2*(1 - normcdf(z))) ;
   
   TABLE1(1,:) = {'VAR','COEFF','Z-VAL','PROB'};

   TABLE1(2:1+length(b),:) = C ;
   
else %otherwise, just report the z-values
   
   TABLE1(1,:) = {'VAR','COEFF','Z-VAL'};
   
   TABLE1(2:1+length(b),:) = C ;
   
end

FINAL_REGRESSION_RESULTS = TABLE1 %Display Regression Results


%Make TABLE2

rho_z = abs(rho)/se(k+3);

TABLE2(1,1) = {'Rho    '};

TABLE2(1,2) = num2cell(rho) ;

TABLE2(1,3) = num2cell(rho_z) ;

if exist('normcdf') == 2 %Check for Stat Toolbox   

   TABLE2(1,4) = num2cell(2*(1 - normcdf(rho_z)));
   
end

AUTOCORRELATION = TABLE2

   
%Make TABLE3

yy = X*b ;

y1 = y - mean(y) ;

y2 = yy - mean(yy) ;

corr = (y1'*y2)/(norm(y1)*norm(y2)) ;

U = ones(n,1) ;

D = eye(n) - (1/n)*U*U' ; %Deviation Matrix

Ps_R2 = (y'*B'*B*X*inv(X'*B'*B*X)*X'*B'*D*B*X*...
        inv(X'*B'*B*X)*X'*B'*B*y)/(y'*B'*D*B*y) ;

TABLE3(:,1) = {'Ps_R2','Sq-Corr','Lik'} ;              
            
TABLE3(:,2) = num2cell([Ps_R2,corr^2,L]) ;
            
GOODNESS_OF_FIT = TABLE3



%Make TABLE4

if exist('chi2cdf') == 2 %Check for Stat Toolbox   
    
   TABLE4(1,:) = {'TEST','VAL','PROB'} ;
   
   TABLE4(2:3,1) = {'LR','Com-LR'};
   
   TABLE4(2:3,2) = num2cell([LR,LR0]);
   
   TABLE4(2,3) = num2cell(1 - chi2cdf(LR,1));
   
   TABLE4(3,3) = num2cell(1 - chi2cdf(LR0,k));
   
else
   
   TABLE4(1,:) = {'TEST','VAL'} ;
   
   TABLE4(2:3,1) = {'LR','Com-LR'};
   
   TABLE4(2:3,2) = num2cell([LR,LR0]);
   
end

TESTS_OF_SAR_MODEL = TABLE4


%Prepare OUTPUT

OUT(:,1) = r0 ; %OLS residuals

OUT(:,2) = W*r0 ; %Weighted residuals

r = B*y - B*X*b ; %SAR Residuals

OUT(:,3) = r ;

OUT(:,4) = W*r ; %Weighted residuals