program define variogram
version 7.0

*	args y

	syntax varname(numeric) [ , INCLude(real 1.5) BWidth(real 0.4) DISCrete SIGMA2(string) RHOOFU(string) ]

	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'" == "")
	if (`include' < 1.0) {
		display " include must be >= 1.0" }
	error (`include' < 1.0)

	preserve

	* Get variance of y 
	qui sum `y'
	local vary = r(Var)

	* Construct pairwise differences and make long
	xtdiff `tis' ulag
	xtdiff `y'   ydif
	by `iis' : gen withinid = _n
	*list `iis' withinid `y' in 1/30
	quietly reshape long ulag ydif , i(`iis' withinid) j(lagf)
	qui drop if (ulag+ydif == .)

	*des
	*list `iis' withinid `y'   lagf ydif in 1/20
	*list `iis' withinid time  lagf ulag   in 1/20

	* Get smooth of pairwise squared differences
	sort ulag
	gen v = 0.5*ydif*ydif

	if ("`discrete'"=="discrete") {
		dis "Computing ANOVA model for v in ulag"
		qui tab ulag , gen(ilag)
		qui regress v ilag* , noconstant
		qui predict vsmth 
		* sum vsmth
	} 
	else {
		dis "Computing smooth lowess model for v in ulag"
		ksmapprox v ulag , num(200) bw(`bwidth') lowess gen(vsmth) nograph
		* sum vsmth
	}

	* Make picture
	qui gen vtmp = v if (v <= `include'*`vary')
	qui summ vtmp
	qui gen vary=`vary'
	scalar `excl' = r(N)
	local exclude = round(100*(_N - `excl')/_N,1) 
	if ("`sigma2'"=="" | "`rhoofu'"=="") {
		# delimit ;
		graph vtmp vsmth vary ulag , c(.ll) s(.ii) 
			ylab xlab 
		t1("Variogram of `y' (`exclude' percent of v_ijk's excluded)")
			l1("v_ijk") b2("| u_ijk |")
			saving(vgram,replace) ;
		# delimit cr
	}
	else {
		qui gen u=ulag
		qui gen sig2fit=`sigma2'
		qui gen gammfit=(`sigma2')*(1 - ( `rhoofu' ))
		# delimit ;
		graph vtmp vsmth vary sig2fit gammfit ulag , 
			c(.lll[-]l[-]) s(.iiii) 
			ylab xlab 
		t1("Variogram of `y' (`exclude' percent of v_ijk's excluded)")
			key3(s(i) c(l) "empirical") 
			key4(s(i) c(l[-]) "fitted") 
			l1("v_ijk") b2("| u_ijk |")
			saving(vgram,replace) ;
		# delimit cr

	}
end
exit