function [RES,RHO,MI] = sar_tom_2(y,X,WW,dist,vnames)

%Variant on sar_tom with fixed set of matrices.

%SAR computes Simultaneous AutoRegression of y on X given weights W 
%    using the multivariate FMINS procedure in MATLAB.

%Written by: TONY E. SMITH, 4/10/98 (Modified: 6/13/99)

%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) WW(:,:,i) = (nxn) matrix i=1:q of spatial weights
%              differing by size of distance radius
%        (iv)  dist = q-vector of distance radii
%        (v)  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'.

%OUTPUT: (1) RES = (n+1)x15 matrix with
%                   RES(1,i) = distance radius

%        (2) RHO = 15x3 matrix with
%                   RHO(i,1) = distance radius
%                   RHO(i,2) = rho value given matrix i
%                   RHO(i,3) = P-value for rho given matrix i

%        (3)  MI = 15x3 matrix with
%                   MI(i,1) = distance radius
%                   MI(i,2) = Moran's I given matrix i
%                   MI(i,3) = P-value for I given matrix i

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


t0 = clock;  %start timer

[n,k] = size(X) ;

q = length(dist);

%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

RES = zeros(n+1,q);  %residuals

RHO = zeros(q,3);  %rho values

MI = zeros(q,3);  %Moran's I


% Start Loop

for count = 1:q      
                       
      RES(1,count) = dist(count);
            
      RHO(count,1) = dist(count);
      
      MI(count,1) = dist(count);      
          
      
      W = sparse(WW(:,:,count));

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

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

info.fmt = '%12.6f'; %print format

%Make Coefficient Display

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

if exist('normcdf') == 2  %check to see if statistics toolbox exists
   
   prob = 2*( 1 - normcdf(abs(z)) );   
    
   
   info.cnames = strvcat('COEFF','Z-VAL','PROB');
   info.rnames = strvcat('VAR','const',vnames);
  
   mat(:,1) = b;
   mat(:,2) = z;
   mat(:,3) = prob;
	
     
else %otherwise, just report the z-values
   
   
   info.cnames = strvcat('COEFF','Z-VAL');
   info.rnames = strvcat('VAR','const',vnames);
  
   mat(:,1) = b;
   mat(:,2) = z;
   
     
end


disp(' ');
disp('REGRESSION CASE:');
disp(' ');

DISTANCE_RADIUS = dist(count)

disp(' ');
disp('FINAL REGRESSION RESULTS:');
disp(' ');

mprint_new(mat,info) %Display Regression Results


%Make Autocorrelation Display

clear( 'mat','info');

info.fmt = '%12.6f';

rho_z = rho/se(k+3);

mat(1,1) = rho;
mat(1,2) = rho_z;

  
info.rnames = strvcat(' ','rho');  

if exist('normcdf') == 2 %Check for Stat Toolbox  
   
   info.cnames = strvcat('VAL','Z-VAL','PROB');

   mat(1,3) = 2*(1 - normcdf(abs(rho_z)));
   
   P_rho = mat(1,3);  %define P-value for later use
   
else
   
   info.cnames = strvcat('VAL','Z-VAL');
   
end



disp(' ');
disp('AUTOCORRELATION RESULTS:');
disp(' ');

mprint_new(mat,info);

   
%Goodness-of-Fit Display

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


clear( 'mat','info');

info.fmt = '%12.6f';


mat(1,1) = Ps_R2;
mat(2,1) = corr^2;
mat(3,1) = L ;

info.rnames = strvcat(' ','Ps_R^2  =','Sq-Corr =','Lik     =');  
   
disp(' ');
disp('GOODNESS-OF-FIT RESULTS:');
disp(' ');

mprint_new(mat,info);


%Test Display

clear('mat','info');

info.fmt = '%12.6f';

mat(1,1) = LR;
mat(2,1) = LR0;

info.rnames = strvcat('TEST','LR','Com-LR');

if exist('chi2cdf') == 2 %Check for Stat Toolbox     
        
   mat(1,2) = 1 - chi2cdf(LR,1);
   mat(2,2) = 1 - chi2cdf(LR0,k);
   
   info.cnames = strvcat('VAL','PROB');
   
else
   
   info.cnames = strvcat('VAL');
      
end

disp(' ');
disp('TESTS OF SAR MODEL:');
disp(' ');

mprint_new(mat,info);

clear('mat','info');

%Prepare OUTPUT

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

RES(2:n+1,count) = r ;

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

RHO(count,2) = rho ;
RHO(count,3) = P_rho;

clear I M u z;

I =  (r'*W*r)/(r'*r) ;

u = ones(n,1);

%Calculate Moran Moments

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

s = u'*W*u ; %Should be equal to n for row normalized W

mu = (n/s)*trace(M*W)/(n-(k+1));

var = ((n/s)^2)*((trace(M*W*M*W') + trace(M*W*M*W) + ...
      ((trace(M*W))^2))/((n-(k+1))*(n-(k+1)+2))) - (mu^2) ;
   
sig = sqrt(var);

z = (I - mu)/sig ;

P = 1 - normcdf(z);

MI(count,2) = I; %Moran value

MI(count,3) = P;  %Moran P-value

end %end i-loop


ELAPSED_TIME = etime(clock,t0)