*! version 1.1.1 PS 16May97.   STB-41 gr27
program define autosmoo
version 4.0
local varlist "req ex min(1) max(2)"
local if "opt"
local in "opt"
local weight "aweight" 
#delimit ;
local options "noGraph GEN(string) GENK(str) GENB(str) Repeat(int 1) LOGit
 Symbol(str) Connect(str) PEn(str) KMIN(int 0) KMAX(int 0) *";
#delimit cr
parse "`*'"
if "`gen'"!="" { confirm new var `gen' }
if "`genk'"!="" { confirm new var `genk' }
if "`genb'"!="" { 
	if "`twice'`mean'`logit'"!="" {
		di in red "genb not available with mean, twice or logit"
		exit 198
	}
	confirm new var `genb' 
}
local nf = 12 +4*("`weight'"!="")
memchk float `nf' double 4 int 3 byte 1
if _rc==90 {exit}
if `repeat'>5 {
	local repeat 5
	noisily di in bl "[repeat set to 5]"
}
parse "`varlist'", parse(" ")
local y `1'
local x `2'
tempvar Y X smooth mse0 mse1 mse2 mse3 sy rsy touse kr kl 
tempvar sxy rsxy beta sx sxx rsx rsxx  kvec
quietly {
	if "`weight'" != "" { 
		tempvar sw w wt
		gen  `w' `exp'
		replace `w'=. if `w'<=0
		local exp "[`weight'=`w']"
	}
	mark `touse' `if' `in'
	markout `touse' `y' `x' `w'
	count if `touse'
	local cnt = _result(1)
	local IN "in 1/`cnt'"
	if `cnt'<5 { noisily error 2001 }
	local med = int(`cnt'/2+.5)
	sum `y' `exp' if `touse'
	local ycen=_result(3)
	gen  `Y'=`y'-`ycen' if `touse'
	_crcslbl `Y' `y'
	if "`x'"!=""{
		gen  `X'=`x' if `touse'
		_crcslbl `X' `x'
	}
	else {
		gen  long `X'=_n if `touse'
		lab var `X' "n"
	}
set graph off
graph `y' `X' `exp' in 1/3, `options' 
set graph on
	sort `X'
	if "`weight'" != "" { 
		gen double `sw'=.
		replace `w'=. if !`touse'
		sum `w'
		noi di in gr "(sum of wgt is " _result(1)*_result(3) ")"
		replace `w' = `w'/_result(3)
		gen `wt' = `w'
		local j 0
		local ks0 "(`sw'[_n+`kr']-cond(_n>`kl',`sw'[_n-`kl'],0)-`wt')"
		local ks "(`sw'[_n+`kr']-cond(_n>`kl',`sw'[_n-`kl'],0))"
		local d3 "`sw'[_N]"
		local wX "`w'*"
		local wtX "`wt'*"
		local wn "`wt'`nk'*"
		local g_swby "by `X':replace `sw'=sum(`w') "   
		local g_wt   "by `X':replace `wt'=`sw'[_N]/_N  "
		local gen_swt "replace `sw'=sum(`wt') `IN'"
	}
	else {
		local ks "(`kr'+`kl')"
		local ks0 "(`kr'+`kl'-1)"
		local d3 "_N"
	}
	gen double `sy'=sum(cond(_n<`cnt',(`X'==`X'[_n-1] |`X'==`X'[_n+1]), /*
		*/           `X'==`X'[_n-1])) `IN'
	local ties=`sy'[`cnt']
	gen `rsy'=sum(`X'!=`X'[_n-1]) `IN'
	local uniq=`rsy'[`cnt']
	if `kmax'== 0 {
		local k3 = int(min(2*`cnt'^.8,`med')+.5)
	}
	else { 	local k3 `kmax' }
	local lim = min(`cnt',100)
	tab `y' in 1/`lim'
	local ny = _result(2)
	if `kmin'== 0 {
		local k0 = int(max(`k3'/250,`cnt'^.2)+.5+`cnt'/(10*`uniq'))
		local k0 = `k0' + int((`cnt'^.2)*((`lim'/`ny')^.2-1)) /*
		*/ + max(0,5-`ny')
	}
	else { 	local k0 `kmin' }
	if `ny'<5 {
		local ny = 1 - `ny'/(4*sqrt(`cnt'))
	}
	else {local ny .95}
	if `k0'>`med' & `k3'>`med' {
		di in red "Can't have kmin greater than `med' (span-min > 1)"
		exit
	}
	local fac = (`k3'/`k0')^(1/3)
	local k1 = int(.5+`k0'*`fac')
	local k2 = int(.5+`k1'*`fac')
	local km = `k0'+1
*	local span=(2*`knn'+1)/`cnt'
	if `ties'> 0 {
	   noi  di "[For span selection only, ties will be randomly separated]"
		tempvar noise
		gen `noise'=`X'-`X'[_n-1] in 2/`cnt'
		sum `noise' if `noise'>0
		replace `noise'=uniform()*_result(5)
	}
	local xcen=(`X'[1]+`X'[`cnt'])/2
	replace `X'=`X'-`xcen' `IN'
	sort `X' `noise'
	if `ties'> 0 { drop `noise' }
        `gen_swt'
	gen  double `sx'=sum(`wtX'`X') `IN'
	gen double `sxx'=sum(`wtX'((`X')^2)) `IN'
	replace `sy'=sum(`wtX'`Y') `IN'
	gen double `sxy'=sum(`wtX'`X'*`Y') `IN'
	gen int `kl'=.
	gen int `kr'=.
	local c0 = `cnt'-3
	local j 3
	while `j' >= 0 {
		local J = `j'+1
		replace `kl' = min(`k`j''+1,_n) `IN'
		replace `kr' = min(`k`j'',`cnt'-_n) `IN'
		diff0 `rsy' `sy' `kr' `kl' `ks0' "`IN'" replace
		if `j'==3 {
		  diff0 `rsx' `sx' `kr' `kl' `ks0' "`IN'" gen
		  diff0 `rsxx' `sxx' `kr' `kl' `ks0' "`IN'" gen
		  diff0 `rsxy' `sxy' `kr' `kl' `ks0' "`IN'" gen
		  local eps = `sxx'[`k2']/(`k2'*1e8)
		  if "`wtX'"!="" {local eps = `sxx'[`k2']/(`sw'[`k2']*1e7) }
		  gen `beta'=(`rsxy'-`rsx'*`rsy')/(`rsxx'-(`rsx')^2) `IN'
		  gen `smooth'=cond((`rsxx'-`rsx'*`rsx') > `eps' , /*
			*/  `rsy'+`beta'*(`X'-`rsx'), `rsy') `IN'
		}
		else if `k`j'' <= `med' & `k`j''<`k`J''{ 
		  local skip
		  replace `smooth' = `rsy' `IN'
		  local lim = `k`j''
		  local N  "_n+`k`j''"
		  if "`wt'"!="" {
			local J  "(`sw'[`N']-`wt')"
		  }
		  else { local J  "(`N' - 1)" }
		  local RSX "(`sx'[`N']-`wtX'`X')"
		  local RSXX "(`sxx'[`N']-`wtX'(`X')^2)"
		  local RSXY "(`sxy'[`N']-`wtX'`X'*`Y')"
		  local RSY "(`sy'[`N']-`wtX'`Y')"
		  replace `beta'= (`RSXY'-`RSX'*`RSY'/`J') /*
		  */	/(`RSXX'-`RSX'^2/`J') in 1/`lim'
		  replace `smooth' = `beta' * (`X'-`RSX'/`J') + `RSY'/`J' /*
		  */	if `beta'!=. in 1/`lim' 
*
		  local lim = `cnt'-`k`j''+1 
		  local N  "_n-`k`j''-1" 
		  if "`wt'"!="" {
			local J "(`sw'[`cnt']-`sw'[`N']-`wt')" 
		  }
		  else { local J  "(_N-_n+`k`j'')" }
		  local RSX "(`sx'[`cnt']-`sx'[`N']-`wtX'`X')"
		  local RSXX "(`sxx'[`cnt']-`sxx'[`N']-`wtX'(`X')^2)"
		  local RSXY "(`sxy'[`cnt']-`sxy'[`N']-`wtX'`X'*`Y')"
		  local RSY "(`sy'[`cnt']-`sy'[`N']-`wtX'`Y')"
		  replace `beta'= (`RSXY'-`RSX'*`RSY'/`J') /*
		  */	/(`RSXX'-`RSX'^2/`J') in `lim'/`cnt'
		  replace `smooth' = `beta' * (`X'-`RSX'/`J') + `RSY'/`J' /*
		  */   if `beta'!=. in `lim'/`cnt' 
		}
		else {local skip "S"}
*noi gr `smooth' `X'
		if "`skip'"!="" {
		  gen `mse`j''=.
		} 
		else {
		  local B0 = (`y'[`med']-`ycen'-`smooth'[`med'])^2
		  replace `beta'=sum((`y'-`ycen'-`smooth')^2 - `B0') `IN'
		  gen `mse`j''= `beta'[min(`cnt',_n+`km')] /*
		  */	   - cond(_n<`km'+1,0,`beta'[_n-`km'-1]) /*
		  */	+ (min(`cnt',_n+`km')-max(0,_n-`km'-1))*`B0' `IN'
		}
		local j = `j' -1
	}
	replace `beta'=min(`mse0',`mse1',`mse2',`mse3')/(`ny') `IN'
	gen int `kvec'= cond(`mse3'<=`beta',`k3',/* 
	*/	cond(`mse2'<=`beta',`k2',/* 
	*/	cond(`mse1'<=`beta',`k1',`k0' ))) `IN'
	replace `beta'=sum(log(`kvec')) `IN'
	replace `kl'=min(int(`cnt'^.8/4+5)+1,_n) `IN'
	replace `kr'=min(int(`cnt'^.8/4+5),`cnt'-_n) `IN'
	replace `kvec'=int(exp((`beta'[_n+`kr']-  /*
	*/	cond(_n>`kl',`beta'[_n-`kl'],0))/(`kl'+`kr'))+.5) `IN'
	if `ties'>0 {
		`g_swby'	/* sw must be defined before d3 is used */
		`g_wt'
		by `X':replace `sy'=sum(`wX'`Y') if `touse'  
		by `X':replace `smooth'=`sy'[_N]/`d3' if `touse' 
		by `X':replace `kvec'=sum(`kvec') if `touse' 
		by `X':replace `kvec'=int(.5+`kvec'[_N]/_N) if `touse' 
	}
	else {	replace `smooth'=`Y'	 `IN' }
	`gen_swt'
	local kk  int(`kvec'/sqrt(`repeat') +.5)
	replace `kl' = min(`kk'+1,_n) `IN'
	replace `kr' = min(`kk',`cnt'-_n) `IN'
	diff `rsx' `sx' `kr' `kl' `ks' "`IN'" replace
	diff `rsxx' `sxx' `kr' `kl' `ks' "`IN'" replace
	while `repeat'>0 {
		replace `sy'=sum(`wtX'`smooth') `IN'
		diff `rsy' `sy' `kr' `kl' `ks' "`IN'" replace
		replace `sxy'=sum(`wtX'`X'*`smooth') `IN'
		diff `rsxy' `sxy' `kr' `kl' `ks' "`IN'" replace
		replace `beta'=(`rsxy'-`rsx'*`rsy')/(`rsxx'-(`rsx')^2) `IN'
		replace `smooth'=cond((`rsxx'-`rsx'*`rsx') > max(0, /*
			*/ `rsxx'/1e8), `rsy'+`beta'*(`X'-`rsx'), `rsy') `IN'
		local repeat=`repeat'-1
	}
	if `ties'>0 {
		by `X':replace `sy'=sum(`smooth')  if `touse'
		by `X':replace `smooth'=`sy'[_N]/_N   if `touse'
	}
	replace `X'=`X'+`xcen'	 `IN'
	replace `Y'=`Y'+`ycen'	 `IN'
	replace `smooth'=`smooth'+`ycen' `IN'	
	if "`logit'"!="" {
		sum `y' if `touse'
		if _result(5)<-5e-9 | _result(6)>1+5e-9 {
			nois di "[Yvar out of range [0,1]. logit not available]"
		} 
		else {
		sum `kl'
		local mkl=2*_result(6)
		sum `y' if `y'>0 & `y'<1 &`touse'
		local aL = min(1/(min(`mkl',1)*`cnt'+1),_result(5))
		local aU = max(1-1/(min(`mkl',1)*`cnt'+1),_result(6))
		replace `smooth'=log(cond(`smooth'< `aL', `aL'/(1-`aL'), /*
		 */  cond(`smooth'>`aU',`aU'/(1-`aU'),`smooth'/(1-`smooth')) /*
		 */	))  if `smooth' !=. `IN'
		local laL = log(`aL'/(1-`aL'))
		local laU = log(`aU'/(1-`aU'))
		local adj=(`laU'-`laL')/25
		replace `Y'=cond(`y'<`aL', `laL'-(1.3+uniform())*`adj', /*
		 */  	    cond(`y'>`aU', `laU'+(1.3+uniform())*`adj', /*
		 */    		log(`y'/(1-`y')) )) `IN'
		local y `Y'
		}
		drop `sx' `sxx' `sxy' `sy'
	}
* GRAPH (If required)
	if "`graph'"!="nograph" {
		if "`symbol'" ==""{ 
			if `cnt'<300 {	local symbol o }
			else         {	local symbol . }
		}
		local symbol  "`symbol'i" 
		if "`connect'"==""{ local connect ".l" }
		if "`pen'"    ==""{ local pen     "23`pp'" }
		noisily graph `y' `smooth' `l' `u' `X' `exp', `options' /*
		 */ s(`symbol') c(`connect') pen(`pen')
	}
	if "`gen'"  !="" {rename `smooth' `gen'}
	if "`genk'" !="" {rename `kvec' `genk'}
	if "`genb'" !="" {rename `beta' `genb'}
}
end

prog def diff0 
	local new `1'
	local old `2'
	local kr `3'
	local kl `4'
	local ks `5'
	local IN `6'
	local com `7'
   	`com' `new' =(`old'[_n+`kr'] - cond(_n>`kl',`old'[_n-`kl'],0)   /*
	*/	- `old' + cond(_n>1,`old'[_n-1],0)	)/`ks'  `IN'
end 

prog def diff 
	local new `1'
	local old `2'
	local kr `3'
	local kl `4'
	local ks `5'
	local IN `6'
	local com `7'
   	`com' `new' =(`old'[_n+`kr'] - cond(_n>`kl',`old'[_n-`kl'],0))/`ks' `IN'
end 

*! version 1.0  17 Jun 1997
program define memchk
	version 5.0
	quietly describe, detail short 
	local width = _result(3)
	local ws = _result(6)
	while "`1'"!="" { 
		if "`2'"=="" { error 198 }
		confirm integer number `2' 
		if "`1'"=="int" {
			local width = `width' + 2*`2'
		}
		else if "`1'"=="byte" { 
			local width = `width' + `2'
			}
		else if "`1'"=="long" | "`1'"=="float" {
			local width = `width' + 4*`2'
			}
		else if "`1'"=="double" { 
			local width = `width' + 8*`2'
			}
		else if substr("`1'",1,3)=="str" { 
			local len = substr("`1'",4,.)
			confirm integer number `len'
			local width = `width' + `len'*`2'
		}
		else { 	error 198 } 
		mac shift 2
	}
	if `width' > `ws' { 
		di _col(8) in red "insufficient memory" 
		exit 900
	}

end

exit