function DAT = sac_perm(z,W,n)

%SAC_PERM performs a random-permutation test of Spatial 
%Autocorrelation (SAC) for a given set of residuals
%using Moran's I as well as the least-squares estimate 
%of rho, the product-moment correlation statistic, r,
%and the spatial Durbin-Watson statistic, d.


%Written by: TONY E. SMITH, 5/19/00

%INPUTS: (i)    z = Nx1 vector of data values for test
%        (ii)   W = NxN weight matrix of spatial relations
%                   (assumed to be ROW STANDARDIZED)
%        (iii)  n = number of simulations to be run
%
%OUTPUT: DAT = (n+1,4)-matrix of index values. 
%              (i)  First row contains actual values based on z.
%              (ii) Last n rows contain the random-sample values 

%SCREEN OUTPUTS: (1) Table showing the range (MIN,MAX) of random sample 
%                    values for each of the four indices (I, rho, r, D-W). 
%                (2) Table of actual (I, rho, r, D-W) values for data
%                    listing significance (P-value) for each index.
%                    Low values of "SIG" indicate the presence of 
%                    significant spatial autocorrelation.

%NOTES: (i) For sake of generality, data z is transformed to deviation form. 
%       This is not necessary for OLS residuals (which already sum to zero)
%       But is necessary for other types of data.
%        
%       (ii) If ties occur, .5 is added to the order counts, thus placing the 
%       ranking of the observed value in the middle of all tied values  


N = length(z);

u = ones(N,1);

z = z - mean(z)*u;  %Put z in deviation form 

%Compute Initial Stats

I0 = (z'*W*z)/(z'*z) ;

rho_0 = (z'*W*z)/(z'*W'*W*z) ;

r_0 = (z'*W*z)/(norm(z)*norm(W*z)) ;

d_0 = (norm(z - W*z)/norm(z))^2 ;

Cnt = zeros(4,1); %vector of counts

v = [I0,rho_0,r_0,d_0]' ;

DAT(1,:) = v';


%Compute test values

h = waitbar(0,'Start Simulations...');

for i = 1:n
   
   %Compute random permutation of locations
   
   P = randperm(N); %MatLab Routine
   
   %Compute Stats
   
   I = (z(P)'*W*z(P))/(z(P)'*z(P)) ;
   
   rho = (z(P)'*W*z(P))/(z(P)'*W'*W*z(P)) ;

   r = (z(P)'*W*z(P))/(norm(z(P))*norm(W*z(P))) ;

   d = (norm(z(P) - W*z(P))/norm(z(P)))^2 ;
   
   DAT(i+1,:) = [I,rho,r,d];
   
   if I > I0  %compare Moran values
       
       Cnt(1) = Cnt(1) + 1;
       
   elseif I == I0
           
       Cnt(1) = Cnt(1) + .5;
       
   end
    
   if rho > rho_0  %compare rho values
       
       Cnt(2) = Cnt(2) + 1;
       
   elseif rho == rho_0
           
       Cnt(2) = Cnt(2) + .5;
       
   end
           
   if r > r_0  %compare correlation values
       
       Cnt(3) = Cnt(3) + 1;
       
   elseif r == r_0
           
       Cnt(3) = Cnt(3) + .5;
       
   end
   
   if d < d_0  %compare Durbin-Watson values
       
       Cnt(4) = Cnt(4) + 1;
       
   elseif d == d_0
           
       Cnt(4) = Cnt(4) + .5;
       
   end   
   
   waitbar(i/n)
   
end %end for

close(h)

%Construct P-value

P_vals = (Cnt+1)/(n+1); %add one to Cnt [so Cnt(i)=0 yields Pi = 1/(n+1)]

C = zeros(4,2);

C(:,1) = DAT(1,:)';

C(:,2) = P_vals;
   
%MAKE SCREEN OUTPUT

M = DAT(2:n+1,:);

MAX = max(M); 
   
MIN = min(M);  
   
%Make Value-Range Display
   
info.cnames = strvcat('Moran','rho','corr','D-W');
info.rnames = strvcat('INDEX','MAX','MIN');
  
mat(1,:) = MAX;
mat(2,:) = MIN;

disp(' ');
disp('RANGE OF RANDOM-PERMUTATION INDEX VALUES:');
disp(' ');

mprint_new(mat,info) %Display Regression Results

%Make Table of Results

info.rnames = strvcat('INDEX','Moran','rho','corr','D-W');
info.cnames = strvcat('VALUE','SIGNIF');
  
disp(' ');
disp('TABLE OF SIGNIFICANCE LEVELS:');
disp(' ');

mprint_new(C,info) %Display Regression Results