function [Pval_1,Pval_2] = moran_sim(X,W,I0,R0,sims)

%MORAN_SIM estimates the actual distribution of Moran's I and Rho under
%the null hypothesis of independence. This is used to calculated actual
%P-values for observed I and Rho and compare their distributions.

%INPUTS:  X = (n x k) matrix of explanatory variables
%         W = (n x n) weight matrix
%        I0 = observed value of I
%        R0 = observed value of Rho
%      sims = number of simulations (default = 10000)
%
%OUTPUT; Pval_1 = P-value for I where by assumption, small P denotes
%               significant positive spatial autocorrelation.
%        Pval_2 = P-value for Rho
%
%SCREEN OUTPUT: Histogram of I-distribution plus a red verticle line
%               to indicate the position of I0.Same for Rho.


if nargin == 4
    
    sims = 10000;
    
end

[n,k] = size(X);

X = [ones(n,1),X]; %add unit vector

M = eye(n) - X*inv(X'*X)*X';

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

I = zeros(sims,1);

Rho = zeros(sims,1);

G = M*W*M;

count_1 = 0;  %count of I-excedences

count_2 = 0; %count of Rho-excedences

for i = 1:sims
    
    e = E(:,i);
    
    num = e'*G*e;
    
    I(i) = num/(e'*M*e);
    
    if I(i) >= I0
        
        count_1 = count_1 + 1;
        
    end
    
    Rho(i) = num/(e'*M*W'*W*M*e);
    
    if Rho(i) >= I0
        
        count_2 = count_2 + 1;
        
    end
    
end

%Calculate Pval

Pval_1 = (count_1 +1)/(sims + 1);

mu_1= mean(I);

mu = trace(M*W*M)/(n - (k + 1));

disp('***************************************');
disp(' ');
disp(['P-value for I = ',num2str(Pval_1)]);
disp(' ');
disp(['simulated mean value of I = ',num2str(mu_1)]);
disp(' ');
disp(['theoretical mean value of I = ',num2str(mu)]);
disp(' ');

mu_2 = mean(Rho);

Pval_2 = (count_2 +1)/(sims + 1);

disp('***************************************');
disp(' ');
disp(['P-value for Rho = ',num2str(Pval_2)]);
disp(' ');
disp(['simulated mean value of Rho = ',num2str(mu_2)]);
disp(' ');

N = hist_ts(I,25); 

y = [0 max(N)/2];

x = I0*ones(1,2);

hist_ts(I,25); hold on;

plot(x,y,'r','LineWidth',5); 

title('Distribution of I Values','Fontsize',12, 'FontWeight', 'bold');

pause;

clf;

hist_ts(Rho,25); hold on;

x = R0*ones(1,2);

plot(x,y,'r','LineWidth',5); 

title('Distribution of Rho Values','Fontsize',12, 'FontWeight', 'bold');