function [I,P] = moran(y,X,W)

%MORAN computes the classical Moran statistic and uses the 
%normal approxation to determine significance.

%Written by: TONY E. SMITH, 2/27/00

%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 (must NOT be
%                    the identity matrix, or degeneracy occurs)


%OUTPUTS: (i)   I = Moran's I-value
%         (ii)  P = p-value (for positive correlation)
%         (iii) r = n:1 vector of OLS residuals

[n,k] = size(X) ;

%Augment X with unit vector

X0 = X ;

clear('X') ;

u = ones(n,1);

X(:,1) = u ;

X(:,2:k+1) = X0 ;

%Initial OLS estimation of beta

b = X\y ;

r = y - X * b ; %OLS residuals

I = (r'*W*r)/(r'*r); %Moran's I

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