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) structure of additional settings
%                opts.maxdist = (approximate) maximum distance
%                             to be considered. (Default = Dmax/2)
%                opts.binmax = maximum number of bins to be used.
%                            (Default = 100)
%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

dbstop if error

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

e0 = y - XX*b ; %OLS residuals

%Construction of initial variogram estimates

M(:,3) = e0 ;

disp(' ');
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

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
    
   if rank(C) < n
       
       C = C + .0001*eye(n); %Avoid almost singular matrices
       
       disp('C augmented to avoid almost singularity');
       
   end

   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);
   
   disp(' ');
   disp(['del = ',num2str(del)]);
     
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');

%Check for flat variogram

if p(2) - p(3) < 10^(-5)
    
    p(1) = 0; 
    
    disp(' ');
	disp('*************************************************');
	disp(' ');
	disp('VARIOGRAM is essentially FLAT. RANGE set to ZERO');
	disp(' ');
	disp('**************************************************');
    
end

mat(:,1) = p';


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

mprint_new(mat,info);

%Plot final variogram

DAT{1} = p;
DAT{2} = results.dist;
DAT{3} = results.vgm;

var_plot(DAT);

d = opts.maxdist;

disp(['MAXDIST = ',num2str(d)]) ;  %Print MaxDist to Screen
disp(' ');


%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} = e0;
OUT{3} = y - XX*b;
OUT{4} = AA*(y - XX*b) ; %Final Regression Residuals
OUT{5} = NEW_DAT ;