function [RHO,BETA,V] = sar_sim(rho,v,b,X,W,sims,E)

%SAR_SIM simulates a spatial autoregressive process with
%        all parameters and X,W data given
%        
%INPUTS: (i)     rho = autoregressive parameter
%        (ii)      v = variance
%        (iii)     b = beta vector
%        (iv)      X = (n x k) attribute matrix
%        (v)       W = (n x n) weight matrix
%        (vi)   sims = number of simulations
%        (vii)     E = (optional) Eigenvalues of W
%
%OUTPUT: RHO    = vector of rho estimates
%         V     = vector of variance estimates
%         BETA  = (k+1 x sims) matrix of beta estimates.
%         Ps_R2 = Pseudo R-square (adjusted)
%         R2    = R-square (adjusted)

[n,k] = size(X);

%Augment X with unit vector

X0 = X ;

clear('X') ;

X(:,1) = ones(n,1) ;

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

b0 = b;

%Compute inv(I - rho*W)

B0 = eye(n) - rho*W;

B0inv = inv(B0);

sig = sqrt(v);

%Compute eigenvalues

val = real(eig(full(W)));
      
[II,JJ,V] = find(val); %find nonzero values

v0 = min(V);

v1 = max(V);

if v0 ~= 0
    
    low = 1/real(v0) ;
    
else
    
    low = 0;
    
end

high = 1/v1 ;

%Simulate data

Z = normrnd(0,1,n,sims) ;

RHO = zeros(sims,1);

V = zeros(sims,1);

BETA = zeros(k+1,sims);

Ps_R2 = zeros(sims,1);

Ps_R2_true = zeros(sims,1);

OLS_R2 = zeros(sims,1);

OLS_R2_true = zeros(sims,1);

R2_new = zeros(sims,1);

U = ones(n,1) ;

D = eye(n) - (1/n)*U*U' ; %Deviation Matrix		

%START SIMULATIONS

for i = 1:sims
   
   y = X*b0 + sig*B0inv*Z(:,i);
   
   %*****************************************************************
   %First estimate residuals under beta = 0 model for goodness of fit   
   %******************************************************************
   
    criterion = .0001; %internal convergence critirion
	itermax = 100;     %internal iteration limit
	converge = 1;
	
    e0 = y - mean(y);    
	econverge = e0;
	iter = 1;
    	
    u = ones(n,1);
    
	while (converge > criterion & iter < itermax) 
       
       rho = fminbnd('sar_lik',low,high,[],low,high,val,W,e0,u) ;
       
		B = speye(n) - rho * sparse(W) ;
	
		b = inv(u'*B'*B*u)*u'*B'*B*y;
        
              
%         bb = (mean(y) - (rho/n)*u'*W*y)/(1-rho);
	
		e = y - u * b ;  
       
        converge = sum(abs(e - econverge));
	
		econverge = e ;
	
		iter = iter + 1;
	
	end %end estimation loop
	
	if (iter == itermax)
       
       error('No convergence within 100 iterations');
       
	end
   
    e0 = e;	%Null Residuals	
    
   %*********************************
   %Now estimate SAR Model
   %*********************************
   
   %Initial OLS estimation of beta
   
   b = X\y ;
   
   r0 = y - X*b ; %OLS residuals
   
   var0 = (1/n)*r0'*r0 ; %OLS error variance   
      
   %Estimate rho and beta
   
    criterion = .0001; %internal convergence critirion
	itermax = 100;     %internal iteration limit
	converge = 1;
	
	econverge = r0;
	iter = 1;
	e = r0;
	
	while (converge > criterion & iter < itermax) 
       
       rho = fminbnd('sar_lik',low,high,[],low,high,val,W,e,X) ;
       
		B = speye(n) - rho * sparse(W) ;
	
		b = inv(X'*B'*B*X)*X'*B'*B*y;
	
		e = y - X * b ;  
       
        converge = sum(abs(e - econverge));
	
		econverge = e ;
	
		iter = iter + 1;
	
	end %end estimation loop
	
	if (iter == itermax)
       
       error('No convergence within 100 iterations');
       
	end
   
   RHO(i) = rho;
   
   V(i) = (1/n) * e'*B'*B*e;
   
   BETA(:,i) = b;   
   
      
	%Compute PSEUDO R-SQUARE (Uses Transformed Model)
	
	e = B*y - B*X*b;	%Transformed Residuals
	
	Ps_R2(i) = 1 - ((n-1)/(n-(k+1)))*(e'*D*e)/(y'*B'*D*B*y);
    
    e = B0*y - B0*X*b0;  %True Transformed Residuals
    
    Ps_R2_true(i) = 1 - ((n-1)/(n-(k+1)))*(e'*e)/(y'*B0'*D*B0*y);
    
    e = y - X*b;	%Raw Residuals
    
    R2_new(i) =  1 - ((n-1)/(n-(k+1)))*(e'*e)/(e0'*e0);	
      
	OLS_R2(i) = 1 - ((n-1)/(n-(k+1)))*(e'*D*e)/(y'*D*y);  
    
    e = y - X*b0;	%True Raw Residuals
    
    OLS_R2_true(i) = 1 - ((n-1)/(n-(k+1)))*(e'*D*e)/(y'*D*y); 
       
end

rho = mean(RHO);

disp(['rho_mean = ',num2str(rho)]);

v = mean(V);

disp(['v_mean = ',num2str(v)]);


b = mean(BETA')