program define acf
version 9.0

** To consider: Should there be an option to remove the mean from each
**              discrete time?  The reason this could be important is
**              that the mean at a given lag of each pair could
**              involve data collected at different times.  Should it
**              be the user's responsibilty to remove trends first?

**              Should the function also report means and standard
**              deviations at each time?

*	args y

	syntax varname(numeric) [ , ]

	tempname excl

	* Get xt variables
	local y `varlist'
	local iis : char _dta[iis]
	local tis : char _dta[tis]

	* Check for correctness
	if ("`iis'" == "" | "`tis'" == "") { 
		display " Wrong! You need to set iis and tis"
              }
	error ("`iis'" == ""  | "`tis'" == "")

	preserve

	* Get variance of y 
	qui xtsumcorr `y'
	local vary = r(var_oc)

	* Construct pairs of variables as fn of lag and make long
	xtdiff `tis' ulag
	xtdiff `y'   ydif
        xtprsum `y'  ysum
	by `iis' : gen withinid = _n
	*list `iis' withinid `y' in 1/30
	quietly reshape long ulag ydif ysum , i(`iis' withinid) j(lagf)
	qui drop if (ulag+ydif == .)
        qui gen yt  =(ysum-ydif)/2   /* y at a given time t */
        qui gen ytpl=(ysum+ydif)/2   /* y at time t plus lag l */
        * list in 1/4

        * Get pairwise correlations
        qui gen lllag=.
        qui gen lln=.
        qui gen llcor=.
        sort ulag
        scalar cnt=0
        qui levelsof ulag , local(alllagvals)
        foreach lag of numlist `alllagvals' {
          scalar cnt=cnt+1
          qui cor yt ytpl if (ulag==`lag')
          qui replace lllag=`lag' if (_n==cnt)
          qui replace lln=r(N) if (_n==cnt)
          qui replace llcor=r(rho) if (_n==cnt)
        }
        qui gen lltol=2/sqrt(lln-3)

        display "Autocorrelation function of `y' "
        list lllag lln lltol llcor if (lllag!=.)
        * Make picture
        # delimit ;
        graph7 llcor lltol lllag , c(ll[-]) s(oo)
               key1(s(o) c(l)    "autocor fn")
               key2(s(o) c(l[-]) "95% tol limit")
               ylab xlab
               yscale(0,1)
               title("Autocorrelation function of `y'")
               l1("correlation") b2("lag")
               saving(acfun,replace)
               ;
        # delimit cr

end
exit