function results = sem2(y,x,W,info)
% PURPOSE: computes spatial error model estimates
%           y = XB + u,  u = p*W*u + e, using sparse algorithms
% ---------------------------------------------------
%  USAGE: results = sem(y,x,W,info)
%  where: y = dependent variable vector
%         x = independent variables matrix
%         W = contiguity matrix (standardized)
%       info = an (optional) structure variable with input options:
%       info.rmin  = (optional) minimum value of rho to use in search (default = -0.99) 
%       info.rmax  = (optional) maximum value of rho to use in search (default = +0.99)   
%       info.eigs  = 0 to compute rmin/rmax using eigenvalues, (1 = don't compute default)
%       info.convg = (optional) convergence criterion (default = 1e-4)
%       info.maxit = (optional) maximum # of iterations (default = 500)
%       info.lflag = 0 for full lndet computation (default = 1, fastest)
%                  = 1 for MC lndet approximation (fast for very large problems)
%                  = 2 for Spline lndet approximation (medium speed)
%       info.order = order to use with info.lflag = 1 option (default = 50)
%       info.iter  = iterations to use with info.lflag = 1 option (default = 30)     
%       info.lndet = a matrix returned by sar, sar_g, sarp_g, etc.
%                    containing log-determinant information to save time
% ---------------------------------------------------
%  RETURNS: a structure
%         results.meth  = 'sem'
%         results.beta  = bhat
%         results.rho   = rho (p above)
%         results.tstat = asymp t-stats (last entry is rho)
%         results.yhat  = yhat
%         results.resid = residuals
%         results.sige  = sige = e'(I-p*W)'*(I-p*W)*e/n
%         results.rsqr  = rsquared
%         results.rbar  = rbar-squared
%         results.lik   = log likelihood
%         results.nobs  = nobs
%         results.nvar  = nvars (includes lam)
%         results.y     = y data vector
%         results.iter  = # of iterations taken
%         results.rmax  = 1/max eigenvalue of W (or rmax if input)
%         results.rmin  = 1/min eigenvalue of W (or rmin if input)
%         results.lflag = lflag from input
%         results.liter = info.iter option from input
%         results.order = info.order option from input
%         results.limit = matrix of [rho lower95,logdet approx, upper95] intervals
%                         for the case of lflag = 1
%         results.time1 = time for log determinant calcluation
%         results.time2 = time for eigenvalue calculation
%         results.time3 = time for hessian or information matrix calculation
%         results.time4 = time for optimization
%         results.time  = total time taken         
%         results.lndet = a matrix containing log-determinant information
%                          (for use in later function calls to save time)
%  --------------------------------------------------
%  NOTES: if you use lflag = 1 or 2, info.rmin will be set = -1 
%                                    info.rmax will be set = 1
%         For n < 1000 you should use lflag = 0 to get exact results
% --------------------------------------------------    
%  SEE ALSO: prt(results), sar, sdm, sac, far
% ---------------------------------------------------
% REFERENCES: LeSage and Pace (2009), pp. 50-51
% Introduction to Spatial Econometrics (1988)
% For lndet information see: Ronald Barry and R. Kelley Pace, 
% "A Monte Carlo Estimator of the Log Determinant of Large Sparse Matrices", 
% Linear Algebra and its Applications", Volume 289, Number 1-3, 1999, pp. 41-54.
% and: R. Kelley Pace and Ronald P. Barry 
% "Simulating Mixed Regressive Spatially autoregressive Estimators", 
% Computational Statistics, 1998, Vol. 13, pp. 397-418.
% ---------------------------------------------------


% originally written by:
% James P. LeSage, 4/2002
% Dept of Finance & Economics
% Texas State University-San Marcos
% 601 University Drive
% San Marcos, TX 78666
% jlesage@spatial-econometrics.com

% revised by Roger Bivand, March, 2010
% to use moments matrix approach from pages 50-51 in
% LeSage and Pace (2009)


time1 = 0; 
time2 = 0;
time3 = 0;
timet = clock; % start the clock for overall timing

n = length(y);

if sum(x(:,1)) ~= n
    tst = sum(x); % we may have no intercept term
    ind = find(tst == n); % we do have an intercept term
        if length(ind) > 0
            error('sem: intercept term must be in first column of the x-matrix');
        elseif length(ind) == 0 % case of no intercept term
            cflag = 0;
            p = size(x,2);
        end;

elseif sum(x(:,1)) == n % we have an intercept in the right place
     cflag = 1;
     p = size(x,2)-1;
end;

results.cflag = cflag;
results.p = p;


% check size of user inputs for comformability
[n nvar] = size(x);
results.meth = 'sem';

[n1 n2] = size(W);

if n1 ~= n2
error('sem: wrong size weight matrix W');
elseif n1 ~= n
error('sem: wrong size weight matrix W');
end;



% return the easy stuff

results.y = y;
results.nobs = n;
results.nvar = nvar; 

% if we have no options, invoke defaults

if nargin == 3
    info.lflag = 1;
end;



% parse input options
[rmin,rmax,convg,maxit,detval,ldetflag,eflag,order,miter,options] = sem_parse(info);


% compute eigenvalues or limits
[rmin,rmax,time2] = sem_eigs(eflag,W,rmin,rmax,n);


results.rmin = rmin;
results.rmax = rmax;
results.lflag = ldetflag;
results.miter = miter;
results.order = order;

% do log-det calculations
[detval,time1] = sem_lndet(ldetflag,W,rmin,rmax,detval,order,miter);


t0 = clock;
Wx = sparse(W)*x;
Wy = sparse(W)*y;
xpx = x'*x;
xpWx = x'*Wx;
xpWpx = Wx'*x;
xpWpWx = Wx'*Wx;
xpy = x'*y;
xpWy = x'*Wy;
xpWpy = Wx'*y;
xpWpWy = Wx'*Wy;
ypy = y'*y;
ypWy = y'*Wy;
ypWpy = Wy'*y;
ypWpWy = Wy'*Wy;

% do vectorized calculations
rgrid = detval(:,1);
nrho = length(rgrid);
epe = zeros(nrho,1);

% modified to do loop /* RSB */

for i=1:nrho;

 rho = rgrid(i,1);

 Axx = xpx - rho*xpWx - rho*xpWpx + rho*rho*xpWpWx;

 Axy = xpy - rho*xpWy - rho*xpWpy + rho*rho*xpWpWy;

 Ayy = ypy - rho*ypWy - rho*ypWpy + rho*rho*ypWpWy;

 b = Axx\Axy;

 epe(i,1) = Ayy - b'*Axx*b;

end;

like =   -detval(:,2) + (n/2)*log(epe);

[tlike,rindex] = min(like);

rho = rgrid(rindex,1);
results.iter = 0;

time4 = etime(clock,t0);

% compute results 
xs = x - rho*Wx;
ys = y - rho*Wy;
b = (xs'*xs)\(xs'*ys);



results.beta= b;
results.rho = rho;
results.yhat = x*results.beta;
e = y - results.yhat;
results.resid = e;
eD = e - rho*sparse(W)*e;
epe = eD'*eD;
results.sige = epe/results.nobs;

sigu = epe;
sige = results.sige;
parm = [results.beta
        results.rho
        results.sige];

results.lik = f2_sem(parm,y,x,W,detval);


if n <= 500, % t-stats using information matrix (Anselin, 1982 pages 

t0 = clock;
B = (speye(n) - rho*sparse(W));
BI = inv(B); WB = W*BI;
pterm = trace(WB'*WB);
xpx = zeros(nvar+2,nvar+2);
xpx(1:nvar,1:nvar) = (1/sige)*x'*B'*B*x;
% rho, rho
xpx(nvar+1,nvar+1) = trace(WB*WB) + pterm;
% sige, sige
xpx(nvar+2,nvar+2) = n/(2*sige*sige);
% rho, sige
xpx(nvar+1,nvar+2) = -(1/sige)*(rho*trace(WB'*WB) - trace(BI'*WB));
xpx(nvar+2,nvar+1) = xpx(nvar+1,nvar+2);
xpxi = invpd(xpx);
tmp = diag(abs(xpxi));
bvec = [results.beta
        results.rho];
results.tstat = bvec./(sqrt(tmp(1:nvar+1,1)));
time3 = etime(clock,t0);

elseif n > 500 

t0 = clock;

hessn = hessian('f2_sem',parm,y,x,W,detval);
xpxi = invpd(-hessn);
xpxi = diag(abs(xpxi(1:nvar+1,1:nvar+1)));
 tmp = [results.beta
        results.rho];

 results.tstat = tmp./sqrt(xpxi);
time3 = etime(clock,t0);

end; % end of t-stat calculations

ym = y - mean(y);
rsqr1 = sigu;
rsqr2 = ym'*ym;
results.rsqr = 1.0 - rsqr1/rsqr2; % r-squared
rsqr1 = rsqr1/(n-nvar);
rsqr2 = rsqr2/(n-1.0);
results.rbar = 1 - (rsqr1/rsqr2); % rbar-squared
results.lndet = detval;

results.time = etime(clock,timet);
results.time1 = time1;
results.time2 = time2;
results.time3 = time3;
results.time4 = time4;


function [rmin,rmax,convg,maxit,detval,ldetflag,eflag,order,iter,options] = sem_parse(info)
% PURPOSE: parses input arguments for far, far_g models
% ---------------------------------------------------
%  USAGE: [rmin,rmax,convg,maxit,detval,ldetflag,eflag,order,iter] =  sem_parse(info)
% where info contains the structure variable with inputs 
% and the outputs are either user-inputs or default values
% ---------------------------------------------------


% set defaults
% if isoctave()
% 
% % /* RSB */
% 
%     options = [0];
% 
% else

    options = optimset('fminbnd');

% end;

options.MaxIter = 500;

eflag = 1;     % default to not computing eigenvalues
ldetflag = 1;  % default to 1999 Pace and Barry MC determinant approx
order = 50;    % there are parameters used by the MC det approx
iter = 30;     % defaults based on Pace and Barry recommendation
rmin = -0.99;     % use -1,1 rho interval as default
rmax = 0.99;
detval = 0;    % just a flag
convg = 0.0001;
maxit = 500;

fields = fieldnames(info);
nf = length(fields);
if nf > 0

 for i=1:nf

    if strcmp(fields{i},'rmin')

        rmin = info.rmin;  eflag = 1;

    elseif strcmp(fields{i},'rmax')

        rmax = info.rmax;  eflag = 1;

    elseif strcmp(fields{i},'eigs')

        eflag = info.eigs; % flag for compute the eigenvalues

    elseif strcmp(fields{i},'convg')

       options.TolFun = info.convg;

    elseif strcmp(fields{i},'maxit')

        options.MaxIter = info.maxit;  

    elseif strcmp(fields{i},'lndet')

    detval = info.lndet;

    ldetflag = -1;

    eflag = 1;

    rmin = detval(1,1);

    nr = length(detval);

    rmax = detval(nr,1);

    elseif strcmp(fields{i},'lflag')

        tst = info.lflag;

        if tst == 0,

        ldetflag = 0; 

        elseif tst == 1,

        ldetflag = 1; 

        elseif tst == 2,

        ldetflag = 2; 

        else

        error('sar: unrecognizable lflag value on input');

        end;

    elseif strcmp(fields{i},'order')

        order = info.order;  

    elseif strcmp(fields{i},'iter')

        iter = info.iter; 

    end;

 end;

 

else, % the user has input a blank info structure

      % so we use the defaults

end; 



function [rmin,rmax,time2] = sem_eigs(eflag,W,rmin,rmax,n);
% PURPOSE: compute the eigenvalues for the weight matrix
% ---------------------------------------------------
%  USAGE: [rmin,rmax,time2] = far_eigs(eflag,W,rmin,rmax,W)
% where eflag is an input flag, W is the weight matrix
%       rmin,rmax may be used as default outputs
% and the outputs are either user-inputs or default values
% ---------------------------------------------------

if eflag == 0
t0 = clock;
opt.tol = 1e-3; opt.disp = 0;
lambda = eigs(sparse(W),speye(n),1,'SR',opt);  
rmin = 1/lambda;   
rmax = 1;
time2 = etime(clock,t0);
else
time2 = 0;
end;


function [detval,time1] = sem_lndet(ldetflag,W,rmin,rmax,detval,order,iter);
% PURPOSE: compute the log determinant |I_n - rho*W|
% using the user-selected (or default) method
% ---------------------------------------------------
%  USAGE: detval = far_lndet(lflag,W,rmin,rmax)
% where eflag,rmin,rmax,W contains input flags 
% and the outputs are either user-inputs or default values
% ---------------------------------------------------


% do lndet approximation calculations if needed

if ldetflag == 0 % no approximation

t0 = clock;    

out = lndetfull(W,rmin,rmax);

time1 = etime(clock,t0);

tt=rmin:.001:rmax; % interpolate a finer grid

outi = interp1(out.rho,out.lndet,tt','spline');

detval = [tt' outi];

    

elseif ldetflag == 1 % use Pace and Barry, 1999 MC approximation



t0 = clock;    

out = lndetmc(order,iter,W,rmin,rmax);

time1 = etime(clock,t0);

results.limit = [out.rho out.lo95 out.lndet out.up95];

tt=rmin:.001:rmax; % interpolate a finer grid

outi = interp1(out.rho,out.lndet,tt','spline');

detval = [tt' outi];



elseif ldetflag == 2 % use Pace and Barry, 1998 spline interpolation



t0 = clock;

out = lndetint(W,rmin,rmax);

time1 = etime(clock,t0);

tt=rmin:.001:rmax; % interpolate a finer grid

outi = interp1(out.rho,out.lndet,tt','spline');

detval = [tt' outi];



elseif ldetflag == -1 % the user fed down a detval matrix

    time1 = 0;

        % check to see if this is right

        if detval == 0

            error('sem: wrong lndet input argument');

        end;

        [n1,n2] = size(detval);

        if n2 ~= 2

            error('sem: wrong sized lndet input argument');

        elseif n1 == 1

            error('sem: wrong sized lndet input argument');

        end;          

end;





function llike = f2_sem(parm,y,x,W,detval);
% PURPOSE: evaluates log-likelihood -- given ML parameters
%  spatial error model using sparse matrix algorithms
% ---------------------------------------------------
%  USAGE:llike = f2_sem(parm,y,X,W,detm)
%  where: parm = vector of maximum likelihood parameters
%                parm(1:k-2,1) = b, parm(k-1,1) = rho, parm(k,1) = sige
%         y    = dependent variable vector (n x 1)
%         X    = explanatory variables matrix (n x k)
%         W    = spatial weight matrix
%         ldet = matrix with [rho log determinant] values
%                computed in sem.m using one of Kelley Pace's routines
% ---------------------------------------------------                                           
%  NOTE: this is really two functions depending
%        on nargin = 3 or nargin = 4 (see the function)
% ---------------------------------------------------
%  RETURNS: a  scalar equal to minus the log-likelihood
%           function value at the ML parameters
%  --------------------------------------------------
%  SEE ALSO: sem, f2_sem2, f_sem
% ---------------------------------------------------



% written by: James P. LeSage 4/2002

% Texas State University-San Marcos

% 601 University Drive

% San Marcos, TX 78666

% jlesage@spatial.econometrics.com



n = length(y);

k = length(parm);

b = parm(1:k-2,1);

rho = parm(k-1,1);

sige = parm(k,1);





gsize = detval(2,1) - detval(1,1);

i1 = find(detval(:,1) <= rho + gsize);

i2 = find(detval(:,1) <= rho - gsize);

i1 = max(i1);

i2 = max(i2);

index = round((i1+i2)/2);

if isempty(index)

index = 1;

end;

detm = detval(index,2);



e = y - x*b;

ed = e - rho*W*e;

epe = ed'*ed;

tmp = 1/(2*sige);



llike =  detm - (n/2)*log(sige) - (n/2)*log(pi) - tmp*epe;



function lik = f_sem(rho,e,W,detval)
% PURPOSE: evaluates concentrated log-likelihood for the 
%  spatial error model using sparse matrix algorithms
% ---------------------------------------------------
%  USAGE:llike = f_sem(lam,eD,W,detm)
%  where: rho  = spatial error parameter
%         e    = begls residuals
%         W    = spatial weight matrix
%         detm =  matrix with [rho log determinant] values
%                computed in sem.m using one of 
%                Kelley Pace's routines           
% ---------------------------------------------------
%  RETURNS: a  scalar equal to minus the log-likelihood
%           function value at the parameter rho
%  --------------------------------------------------
%  NOTE: this is really two functions depending
%        on nargin = 3 or nargin = 4 (see the function)
% ---------------------------------------------------        
%  SEE ALSO: sem, f_far, f_sac, f_sar
% ---------------------------------------------------



% written by: James P. LeSage 1/2000

% Texas State University-San Marcos

% 601 University Drive

% San Marcos, TX 78666

% jlesage@spatial-econometrics.com





gsize = detval(2,1) - detval(1,1);

i1 = find(detval(:,1) <= rho + gsize);

i2 = find(detval(:,1) <= rho - gsize);

i1 = max(i1);

i2 = max(i2);

index = round((i1+i2)/2);

if isempty(index)

index = 1;

end;

detm = detval(index,2);



ed = e - rho*W*e;

n = length(ed);

lik =  (n/2)*log(ed'*ed) - detm;









function H = hessian(f,x,varargin)

% PURPOSE: Computes finite difference Hessian

% -------------------------------------------------------

% Usage:  H = hessian(func,x,varargin)

% Where: func = function name, fval = func(x,varargin)

%           x = vector of parameters (n x 1)

%    varargin = optional arguments passed to the function

% -------------------------------------------------------

% RETURNS:

%           H = finite differnce hessian

% -------------------------------------------------------



% Code from:

% COMPECON toolbox [www4.ncsu.edu/~pfackler]

% documentation modified to fit the format of the Ecoometrics Toolbox

% by James P. LeSage, Dept of Economics

% Texas State University-San Marcos

% 601 University Drive

% San Marcos, TX 78666

% jlesage@spatial-econometrics.com



eps = 1e-5;



n = size(x,1);

fx = feval(f,x,varargin{:});

 

% Compute the stepsize (h)

h = eps.^(1/3)*max(abs(x),1e-2);

xh = x+h;

h = xh-x;    

ee = sparse(1:n,1:n,h,n,n);

 

% Compute forward step 

g = zeros(n,1);

for i=1:n

  g(i) = feval(f,x+ee(:,i),varargin{:});

end

   

H=h*h';

% Compute "double" forward step 

for i=1:n

for j=i:n

  H(i,j) = (feval(f,x+ee(:,i)+ee(:,j),varargin{:})-g(i)-g(j)+fx)/H(i,j);

  H(j,i) = H(i,j);

end

end



