* Program swiped off of F. Dominici Web page -- 14 mar 02

program define autocor

/* This program is designed to create a scatterplot matrix of the
residuals for longitudinal data.  The arguments that the program takes
are the response (y), the covariates (x), a person identifier (id),
and an indicator of time when the response was measured (time).  In
the results window we get the autocorrelation matrix and the
autocorrelation vector under the stationarity assumption.  In the
graph window, we get two plots: the one on the left is the scatterplot
matrix of residuals and the one on the right is the autocorrelation
function.  Each of these plots is saved in c:\autoc.gph and
c:\acf.gph, respectively.
 
syntax is "autocor y time id" */

version 5.0
preserve
set more off
set graphics off

quietly{

sort `3' `2'			/* sort by id and time */ 
regress `1' `2'			/* regress y on time */
predict res, residuals 		/* save residuals */
tabulate `2', gen(temp)	      /** create temp1 = residuals for no lag, temp2 = resids for 1, etc. **/

/* generate rtempi where rtempi is the residual at time i */

for temp*: gen r@=cond(@==1,res,.)    /** get the residual value for no lag only **/ 
for rtemp*: by `3': replace @=sum(@)  /** for each id sum resids at each time for each person **/
for rtemp*: by `3': replace @=@[_N]   /** make each person's resid for each observation at each time **/
for rtemp*: replace @=cond(@==0,.,@)  /** replace with 0 for no lag **/
	
inspect `2'
local maxt=_result(7)		/* maxt is the number of unique time values */

gen uid=cond(`3'~=`3'[_n-1],1,0) /* uid is an indicator of the first obsn
						on a person */

renpfix rtemp time /* time is now the time point variable */

#delimit ;

/* make scatter plot of residuals for each time combination */

graph time1-time`maxt' if uid==1, matrix half s(o) label saving(autoc, replace)
	 t1("Autocorrelation Scatterplot");

#delimit cr


/* calculate pairwise correlation matrix and show how many obsns
 were used to calculate each correlation within the matrix (listed
 below each correlation) */

noisily pwcorr time1-time`maxt' if uid==1

/* assume stationarity: calculate correlation */

stack time1-time`maxt' if uid==1, into(lag)
local r=`maxt'-1
display `r'
display `maxt'
matrix C=J(`r',1,0) /** rX1 matrix of 0's **/
matrix colnames C=acf
matrix T=J(`r',1,1) /** rX1 matrix of 1's **/
matrix colnames T=tau
matrix list T
matrix list C 

local i=1
local maxid=_N/_stack[_N] /** number of obs per time point **/
while `i'<`maxt'	{
	local m=`maxid'*`i' /** contains # of obs times lag being dealt with **/
	gen lag`i'=lag[_n-`m']	/** lag 2 = val in space 2 spaces ago  **/
	corr lag lag`i'
	scalar c`i'=_result(4)
	local i=`i'+1
			
}

for 1-`r', ltype(numeric) nostop: matrix C[@,1]=c@
matrix list C
drop _all

svmat C, names(col)
svmat T, names(col)
replace tau=sum(tau) if acf~=.
replace acf=. if acf==0
}

list acf 				/* acf is autocorrelation function,
					   tau is time lag */

/* plot acf */

graph acf tau, xlab ylab t1("ACF") yscale(0,1) saving(acf, replace)

/* combine 2 plots on same page */

set graphics on

graph using autoc acf

end