/*Filename: 	autocorr.sas														*/
/*Purpose:		Computes Moran's I, Geary's C, and semivariance for specified      */
/*				spatial lags. Then plots them in a correlogram						*/ 
/*Last Update:	FDN 1.18.07															*/ 

options  ls=80 ps=40 ;                                                          

/*start by setting the graphics options --Note that names of windows fonts are given in quotes*/
goptions reset=global gunit=in 
         ftitle="garamond" ftext="garamond" htitle=.3 htext=.2;

*title1 'Maya Correlograms';  ******* Titles for output here;                                                                                 


/* plot the points locations */
proc gplot data=maya; plot sn * we ;run;


/* do a univarite summary of the data field */

proc univariate data=house;
var date_ad_;
histogram Date_ad_ / cfill=grey;
run;

 
/*invoke Proc IML -- Interactive Matrix Language */
proc iml;                                                                       
use maya   ;                   ******** supply dataset name;
read all var{sn we} into c;    ******** names of northing (y) and easting (x) vars;
read all var{Date_ad_} into x;     ******** name of z variable;                       
                                                                                 
nsteps=    10                ; ******** enter number of lags ;                 
                                                                                
                                                                                
lag=   20                    ; ******** enter lag size;                         
                                                                                
n=nrow(x);                                                                      
distmat= j(n,n,0); gc=j(nsteps,1,0); lagd=gc; semivar=gc;   *dimension vars;    
nedges=gc; moran=gc; diffmat=distmat;  cpmat=distmat;                  
meandist=gc;
mean = (sum(x)/n); 
 do j=1 to n by 1;      *get a matrix of distances between data points;         
     do k=1 to n by 1;  *and matrices of squared differences and crossproducts;                   
       distmat[k,j]= ((c[k,]-c[j,])*(c[k,]-c[j,])`)##.5;                        
       diffmat[k,j]= (x[k,]-x[j,])##2; 
	   cpmat[k,j] =  (x[k,]-mean)#(x[j,]-mean); 	 
    end;                                                                        
 end;                                                                           
z = x - ((sum(x))/n);  *deviations from mean;                                   
ss=ssq(z);                                                                      
variance= ss/(n-1);   
print variance; 
                                                                                
 do i=1 to nsteps;      *compute for the i'th lag;  
	lagd[i,]= i#lag;  
    lh= (lagd[i,]-(lag/2));   hh= (lagd[i,]+(lag/2)); print  lh hh;                                     
    w= (distmat > lh) #  (distmat <= hh); * binary weights matrix;                                                                                                              
	sumw= sum(w);
    meandist[i,]=	(sum(w#distmat))/sumw;              *average diatnce between point pairs;                              
    semivar[i,]=	((sum(w#diffmat))/sumw)*.5; 		*semivariance;                                                                              
    gc[i,]= 		semivar[i,]/variance;    			*gearys c;
	moran[i,] = 		((sum(w#cpmat))/sumw)/ (ss/n);     
	nedges[i,]= sumw/2;                                                 
  end;                                                                          
result= semivar||gc|| moran || lagd|| meandist || nedges;                       
label= {semivar gearys_c morans_i lagdist meandist nedges};                       
create autocorr from result [colname=label];                               
append from result ;                                           


data autocorr; set autocorr;
logmeandist = log10(meandist);
logsemivar=log10(semivar);



proc print;                                                                     

symbol1 value=dot h=.2 i=join ;	*specify the symbol to be used in plots;

proc gplot data=autocorr;
			plot gearys_c * meandist=1;                                      
            plot semivar*meandist =1;
			plot morans_i*meandist =1;
			
run;