function [OUT,cov] = sp_lag_wtd(y,X,W,vnames,wts,eigvals)

%This version of SP_LAG reweights data with a vector of weights, wts.

%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 (Modified 6/14/99)
%
%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)  vnames = (optional) 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'.
%        (v)   eigvals (optional) = vector if eigenvalues of W.
%              [NOTE: for large W-matrices it pays to compute
%               eigenvalues ONCE and then use these as inputs]
%
%OUTPUT: OUT = nX8 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 .
%        OUT(:,5) = y vector .
%        OUT(:,6) = yHat vector .
%        OUT(:,7) = Moran's I value
%        OUT(;,8) = P-value for Moran (normal approx)
%
%               I = Inverse of asymptotic cov matrix
%
%        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


%dstop if error

[n,k] = size(X) ;

%Weight Data

if length(wts) ~= n
    
    error('WRONG DIMENSION ON WTS');
    
end

wts = reshape(wts,n,1); %be sure it is a column vector

y = wts.*y;

X0 = X;

X0 = (wts*ones(1,k)).*X0;

X = [wts,X0];

%Initial OLS estimation of beta

b = X\y ;

r0 = y - X * b ;

var0 = (1/n)*r0'*r0 ; %OLS error variance

%Eigen values

if nargin == 6
   
   val = real(eigvals);
   
   v0 = min(val);
   
   v1 = max(val);
   
elseif nargin < 6
   
   val = real(eig(full(W)));

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

	v0 = min(V);

	v1 = max(V) ;

end

del = .00001;

if v0 ~= 0
        
        low = 1/(-abs(v0)) + del;
        
    else 
        
        low = del;
        
    end

high = 1/v1 - del;


%Initial Estimate of lamda

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

e0 = M*y ;

e1 = M*(W*y) ;


z = fminbnd('sp_lag_lik',low,high,[],low,high,val,e0,e1,n) ; 


%Compute Log Likelihood

lam = min([z ,high - .001]) ;

lam = max([lam,low + .001]) ;

W = sparse(W);

B = speye(n) - lam*W;

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

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

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

LogDet = sum(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(B) ;

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

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

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

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

I(k+2,k+2) = n/2 ;

I(k+2,k+3) = 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) = (var^2)*trace(G*(G + G')) + var*(G*X*b)'*(G*X*b) ;

cov = (var^2)*inv(I) ;

se = sqrt(diag(cov)) ;



%SCREEN OUTPUT

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

%Make Coefficient Display

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

%Construct Variable Names

if ~exist('vnames')

    nvar = k;
    
    vnames = [];
    
    for i=1:nvar
        
        tmp = ['Variable ',num2str(i)];
    
        vnames = strvcat(vnames,tmp);
        
    end;
    
 end;

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('FINAL REGRESSION RESULTS:');
disp(' ');

mprint_new(mat,info) %Display Regression Results

disp(['Variance = ',num2str(var)]);

%Make Autocorrelation Display

clear( 'mat','info');

info.fmt = '%12.6f';

lam_z = lam/se(k+3);

mat(1,1) = lam;
mat(1,2) = lam_z;
   
info.rnames = strvcat(' ','lam');  

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

	mat(1,3) = 2*(1 - normcdf(abs(lam_z))); 

else
   
   info.cnames = strvcat('VAL','Z-VAL');
   
end

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

mprint_new(mat,info);

   
%Make Goodness-of-Fit Display

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


%Compute PSEUDO R-SQUARE (Uses Transformed Model)

%    NOTE: Since the reduced model, By = Xb + e , is a standard
%    regression model WITH INTERCEPT, R-square still makes 
%    perfect sense.

e = B*y - X*b;

Ps_R2 = 1 - (e'*e)/(y'*B'*D*B*y); %Note that 1'e = 0 by OLS

Ps_R2_Adj = 1 - ((n-1)/(n-(k+1)))*(e'*e)/(y'*B'*D*B*y);

%NEW R-SQUARE (uses sp_lag beta in standard R^2).

% b = inv(X'*X)*X'*B*y;
% 
% e = B*y - X*b;
% 
% Ps_R2 = 1 - (e'*D*e)/(y'*D*y);
% 
% Ps_R2_Adj = 1 - ((n-1)/(n-(k+1)))*(e'*D*e)/(y'*D*y);

%Compute corrected AIC (AIC_c) [Burnham and Anderson, MODEL SELECTION...,p.66]

K = k + 3; %Number of parameters estimated (b + var +lam)

AIC = 2*(-L + K);

AIC_c = AIC + 2*K*(K+1)/(n-K-1);


clear( 'mat','info');
   
   
disp(' ');
disp('GOODNESS-OF-FIT RESULTS:');
disp(' ');

disp(strvcat(['Pseudo R-Square      = ',num2str(Ps_R2)])); 
disp(strvcat(['Pseudo R-Square Adj  = ',num2str(Ps_R2_Adj)]));
disp(strvcat(['Squared_Correlation  = ',num2str(corr^2)]));
disp(strvcat(['Log Likelihood Value = ',num2str(L)])); 
disp(strvcat(['AIC = ',num2str(AIC)])); 
disp(strvcat(['AIC_corrected = ',num2str(AIC_c)])); 


%Make Test Display

clear('mat','info');

info.fmt = '%12.6f';

mat(1,1) = LR;

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

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

      
end

disp(' ');
disp('TEST OF SP_LAG MODEL:');
disp(' ');

mprint_new(mat,info);


%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

OUT(:,5) = y ;   %dependent variable

OUT(:,6) = inv(B)*X*b ; %mean predicted values


%Calculate Moran Moments

clear M u z;

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

u = ones(n,1);


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

u = ones(n,1);

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 = (Moran - mu)/sig ;

P = 1 - normcdf(z);

OUT(:,7) = Moran; %Moran value

OUT(:,8) = P;  %Moran P-value