function [LAM,BETA,V,MI] = sp_lag_sim(lam,v,b,X,W,sims,E);

%SP_LAG_SIM simulates a spatial lag process with
%        all parameters and X,W data given
%        
%INPUTS: (i)     lam = 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: LAM    = vector of lam estimates
%         V     = vector of variance estimates
%         BETA  = (k+1 x sims) matrix of beta estimates.
%         MI = Moran's I


[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 - lam*W)

B0 = eye(n) - lam*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);
%    
% low = max(-1,1/(-abs(v0))) ;
% 
% high = 1/v1 ;

%Simulate data

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

LAM = zeros(sims,1);

V = zeros(sims,1);

BETA = zeros(k+1,sims);

MI = zeros(sims,1);

Y = zeros(sims,1);

M = eye(n) - X*inv(X'*X)*X'; %Note that ones(n,1) is in X

U = ones(n,1) ;

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

a = n/sum(sum(W));

%START SIMULATIONS

for i = 1:sims
   
   y = B0inv*(X*b0 + Z(:,i));
   
   Y(i) = sum(y);
   
   [b,lam,v] = sp_lag_silent(y,X,W,val);
   
   LAM(i) = lam;
   
   V(i) = v;
   
   BETA(:,i) = b; 
   
   MI(i) = a*(y'*M*W*M*y)/(y'*M*y);
      
end

lam = mean(LAM);

disp(['lam_mean = ',num2str(lam)]);

v = mean(V);

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


b = mean(BETA');

disp(['b_mean = ',num2str(b)]);