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

%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

%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) = SL residual vector, r .
%        OUT(:,4) = Weighted residuals, W*r .
%
%        NOTE: The residual vector, r, is for the reduced
%              OLS model, B*y = X*b + e. Hence, a comparison
%              of the regressions, (r0 on W*r0) and (r on W*r),
%              indicates the degree to which the SL 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-Squared (in reduced model)
%                 Sq-Corr = squared correlation coefficient
%                 Lik = value of log likelihood
%
%         TABLE3: (SPATIAL-LAG PARAMETER)
%                 Lam estimate plus z-value for significance tests
%                 (If MATLAB Statistics Toolbox is present, the P-values
%                  are also displayed as 'PROB').
%
%         TABLE4: (TEST OF SPATIAL-LAG MODEL)
%                 LR  = likelihood ratio for spatial autocorrelation


[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

%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 ;


%Initial Estimate of lamda

M = eye(n) - X*inv(X'*X)*X' ;

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]) ;

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

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

e = y - lam*W*y - X*b ;

var = (1/n) * e'*e ;

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

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

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

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

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

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

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

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

I(1:k+1,k+3) = (1/var)*X'*G*X*b ;

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')) + (1/var)*(G*X*b)'*(G*X*b) ;

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 = 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

lam_z = lam/se(k+3);

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

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

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

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

SPATIAL_LAG = TABLE2



%Make TABLE3

yy = lam*W*y + 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'*X*inv(X'*X)*X'*D*X*inv(X'*X)*X'*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,1) = {'LR'};
   
   TABLE4(2,2) = num2cell(LR);
   
   TABLE4(2,3) = num2cell(1 - chi2cdf(LR,1));   
     
else
   
   TABLE4(1,:) = {'TEST','VAL'} ;
   
   TABLE4(2,1) = {'LR'};
   
   TABLE4(2,2) = num2cell(LR);
   
end

TEST_OF_SPATIAL_LAG_MODEL = TABLE4

%Prepare OUTPUT

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

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

r = y - lam*W*y - X*b ; %SL Residuals

OUT(:,3) = r ;

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