function OUT = geo_regr(y,X,L,vnames,opts)

%GEO_REGR is an iterative geostatistical regression procedure 
%which uses a spherical covariogram to fit residual covariances.

%Written by: TONY E. SMITH, 2/28/98

%Functions called: distance.m, var_spher.m, spher_mat.m
%
%INPUTS: (i)   y = (n,1) vector of dependent variable values
%        (ii)  X = (n,k) data matrix (Xi, i = 1:k) [usually 
%                   includes polynomial functions of locations].
%        (iii) L = (n,2) matrix of locations (Ei,Ni)
%        (iv)  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'.
%        (v)   opts = optional bandwidth for variogram estimation
%
%OUTPUTS: SCREEN: 

%         Scr_1 = (k+1,2) vector of coeffs plus std.errors
%         Scr_2 = (3,1) vector of estimated variogram parameters
%                  [r,s,a]        
%         DATA:
%         OUT{1} = Beta coefficients       
%         OUT{2} = Regression residuals for OLS
%         OUT{3} = Final regression residuals
%         OUT{4} = Regression residuals for transformed problem
%         OUT{5} = transformed (y,X) data which can be used in
%                  a standard (zero-intercept) regression package.
%                  If the original data variables are denoted by x0(=1),
%                  x1,..,xk, and if transformed values are denoted by
%                  by '*', then each row i of OUT{3} has the form
%                  (y*i,x0*i,x1*i,...,xk*i).


%Initial Settings

[n,k] = size(X) ;

U = ones(n,1) ;

XX(:,1) = U ;

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

%Now rescale data to avoid possible instabilities

S = inv(diag(max(abs(XX)))) ;

Z = XX*S ;

%Define Locations and compute distances

M(:,1:2) = L ;

D = distance_mat(L) ;

Dmin = min(min(D));

Dmax = max(max(D));

%Initial OLS estimates and residuals

bb = Z\y ;

b = S*bb ; %This factors out the rescaling effect as seen below:

%EXPLANATION: Z = (XX)*S => bb = inv(Z'*Z)*Z'*y
%                              = inv[S'*(XX)'*(XX)*S]*S'*(XX)'*y
%                              = inv(S)*inv[(XX)'*(XX)]*inv(S')*S'*(XX)*y
%                              = inv(S)*{inv[(XX)'(XX)]*(XX)'*y}
%                              = inv(S)*b
%             => b = S*bb.

e = y - XX*b ;

%Construction of initial variogram estimates

M(:,3) = e ;

disp('Beginning initial covariogram estimate');

%Set BandWidth if not already specified

if (nargin < 5)
   
   if (Dmin < Dmax/2)
   
   opts.maxdist = Dmax/2 ; %maximum relevant distance
   
	else
   
      error('Dmin >= Dmax/2') ; 
      
   end
   
  
end

BANDWIDTH = opts.maxdist;  %Print Bandwidth to Screen


results = var_spher(M,opts) ;

p = results.param ; %parameter vector for variogram

del = 1.0 ;

%Begin Estimation Loop

disp('Beginning main iteration loop');

while del > .001
   
   %Save current beta value
   
   b0 = b ;
   
	%Construction of Covariance Matrix
   
   C = p(2)*U*U' - spher_mat(p,D) + p(3)*eye(n) ;
   
   
	%New b estimates and residuals

   A = chol(C)';
   
   clear('C'); 
   
   AA = A\eye(n) ;
   
   clear('A') ; %Only AA is used from this point on   
   
   bb = (AA*Z)\(AA*y) ;
   
   b = S*bb ; %This factors out the rescaling effect (as above)
        
   M(:,3) = y - XX*b ;
   
   results = var_spher(M,opts) ; 
     
   p = results.param ; 
   
   del = max(abs(b - b0)./b0)    
     
end

%Calculate standard errors

e = y - XX*b ;

s1 = sqrt((1/(n-(k+1)))*e'*AA'*AA*e) ;

s2 = sqrt(diag(S*inv(Z'*AA'*AA*Z)*S)) ; %Again uses rescaling

se = s1*s2 ;


%SCREEN OUTPUT

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

%Make Coefficient Display

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

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

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

mprint_new(mat,info) %Display Regression Results


%Make Variogram Display

clear( 'mat','info');

info.fmt = '%12.6f';

info.rnames = strvcat(' ','Range','Sill','Nugget');

mat(:,1) = p';


disp(' ');
disp('SPHERICAL VARIOGRAM PARAMETERS:');
disp(' ');

mprint_new(mat,info);



%PREPARE DATA OUTPUT

%Summarize Transformed Data

NEW_DAT(:,1) = AA*y ;
NEW_DAT(:,2:k+2) = AA*XX ;

%Prepare DATA OUTPUT

OUT{1} = b;
OUT{2} = e;
OUT{3} = y - XX*b;
OUT{4} = AA*(y - XX*b) ; %Final Regression Residuals
OUT{5} = NEW_DAT ;